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Abstract 


The  focus  of  this  research  is  on  single  platform  geolocation  methods  where  the 
position  of  a  single  stationary  radio  frequency  emitter  is  estimated  from  multiple 
simulated  angle  and  frequency  of  arrival  measurements  taken  from  a  single  moving 
receiver  platform.  The  analysis  scenario  considered  consists  of  a  single  6U  CubeSat  in 
low  earth  orbit  receiving  radio  frequency  signals  from  a  stationary  emitter  located  on  the 
surface  of  the  Earth.  A  multiple  element  receive  antenna  array  and  the  multiple  signal 
classification  algorithm  are  used  to  estimate  the  angles  of  arrival  of  an  impinging  signal. 
The  maximum  likelihood  estimator  is  used  to  estimate  the  frequency  of  arrival  of  the 
received  signal.  Four  geolocation  algorithms  are  developed  and  the  accuracy  performance 
is  compared  to  the  Cramer-Rao  lower  bounds  through  Monte  Carlo  simulations.  Results 
from  a  system  parameter  sensitivity  analysis  show  the  combined  angle  and  frequency  of 
arrival  geolocation  maximum  likelihood  estimator  consistently  outperforms  the  other 
three  geolocation  algorithms. 
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SINGLE  PLATFORM  GEOLOCATION  OF  RADIO  FREQUENCY  EMITTERS 


I.  Introduction 

Geolocation  of  Radio  Frequency  (RF)  emitters  involves  estimating  the  position  of  an 
emitter  from  the  received  RF  signals  at  one  or  more  receivers.  In  the  case  of  passive 
geolocation  [1],  only  the  phenomenology  of  the  received  signals  (angle  of  arrival,  time 
delay,  Doppler  frequency  shift,  etc.)  is  used  to  estimate  the  position  of  the  RF  emitter, 
rather  than  the  emitter  providing  its  position  as  a  message  contained  in  the  RF  signal. 
Geolocation  methods  involving  multiple  coordinated  RF  receivers  include:  Angle  of 
Arrival  (AOA)  [2],  time  difference  of  arrival  [3],  frequency  difference  of  arrival  [1],  and 
direct  position  determination  [4].  Geolocation  methods  using  a  single  moving  receiver 
include  AOA  [3]  and  Frequency  of  Arrival  (FOA)  [5].  The  focus  of  this  thesis  research  is 
the  single  moving  platform  geolocation  methods  available  to  estimate  the  position  of  a 
single  stationary  RF  emitter  using  multiple  angle  and  frequency  measurements  from  the 
received  RF  signals. 

1.1  Problem  Statement  and  Research  Objective 

The  scenario  considered  in  this  research  is  to  geolocate  a  single  stationary  RF  emitter 
using  the  signals  received  by  an  antenna  array  on  a  single  6  Unit  (6U)  CubeSat  [6]  moving 
in  Low  Earth  Orbit  (LEO).  The  goal  of  this  research  is  to  develop,  implement,  analyze, 
and  compare  single  platform  geolocation  algorithms.  A  simulation  framework  is 
developed  to  conduct  a  system  parameter  sensitivity  analysis  to  assess  parameter  impact 
on  the  performance  of  the  geolocation  algorithms. 
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1.2  Underlying  Assumptions 

The  following  underlying  assumptions  apply  to  the  scope  of  the  research.  The  single 
moving  6U  CubeSat  consists  of  a  calibrated  multiple  element  antenna  array,  phase 
coherent  RF  receivers,  and  a  guidance,  navigation,  and  control  subsystem  to  provide 
position,  velocity,  and  attitude  data  to  a  dedicated  payload  processor.  A  6U  CubeSat  is 
defined  [6]  with  the  exterior  dimensions  of  12  x  24  x  36  cm  and  a  total  mass  less  than 
12  kg.  The  single  terrestrial  RF  emitter  is  assumed  to  be  stationary  with  no  co-channel 
interference  from  multiple  emitters.  The  following  topics  are  beyond  the  scope  of  this 
thesis  and  are  left  as  considerations  for  future  research:  implementation  of  specific 
antenna  arrays,  RF  emitters,  and  RF  receivers;  filtering  and  pre-processing  of  the  received 
RF  signals;  multiple  emitter  segregation  through  data  association  of  the  received  RF 
signals;  and  the  effect  of  co-channel  interference  on  geolocation  algorithm  performance. 

1.3  Thesis  Organization 

This  thesis  is  organized  as  follows:  Chapter  II  discusses  existing  angle  of  arrival 
geolocation  algorithms  and  signal  parameter  estimation  methods;  Chapter  III  develops  the 
geolocation  algorithms  and  simulations  used  in  this  thesis;  Chapter  IV  reports  the  results 
of  a  system  parameter  sensitivity  analysis  to  assess  the  impact  on  the  performance  of  the 
various  geolocation  algorithms;  and  Chapter  V  presents  the  overall  conclusions  of  this 
research  with  recommendations  for  future  work. 
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II.  Background 


This  chapter  is  organized  as  follows:  Section  2.1  describes  the  coordinate  systems  and 
reference  frames  used  in  this  thesis;  Section  2.2  details  the  Multiple  Signal 
Classification  (MUSIC)  angle  of  arrival  algorithm  and  associated  theoretical  performance; 
Section  2.3  summarizes  two  existing  angle  of  arrival  geolocation  methods;  Section  2.4 
describes  maximum  likelihood  estimation,  Gauss-Newton  iterations,  and  the  Cramer-Rao 
lower  bound;  Section  2.5  details  the  maximum  likelihood  parameter  estimation  of  a 
sinusoidal  signal;  and  Section  2.6  describes  a  method  to  visualize  the  confidence  of 
estimated  parameters. 

2.1  Coordinate  Systems  and  Reference  Frames 

The  following  coordinate  systems  and  reference  frames  are  used  throughout  this  thesis 
to  define  relative  position  and  attitude  geometries.  The  Earth  Centered  Earth 
Fixed  (ECEF)  coordinate  system  is  used  for  positioning  in  terms  of  3D  Cartesian 
coordinates.  The  Geodetic  coordinate  system  using  the  World  Geodetic  System  (WGS)  84 
model  is  used  for  positioning  in  terms  of  latitude,  longitude,  and  altitude.  The  local  ECEF 
reference  frame  is  a  translated  ECEF  coordinate  system.  The  East  North  Up  (ENU) 
reference  frame  is  used  to  define  the  satellite  local  reference  frame.  The  sensor  reference 
frame  is  used  to  define  the  antenna  array  geometry  with  respect  to  the  satellite  local 
reference  frame.  The  relationship  between  the  ECEF,  geodetic,  and  ENU  coordinate 
systems  is  shown  in  Figure  2.1. 

The  ECEF  frame  is  a  geocentric  right  handed  3D  Cartesian  coordinate  system  with  the 
origin  at  the  center  of  mass  of  the  Earth.  The  Xe-axis  points  towards  the  intersection  of  the 
equator  and  the  prime  meridian  (0°  latitude  0°  longitude),  the  Ze- axis  points  towards  the 
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Figure  2.1:  Relationship  between  the  ECEF  (Xe,Ye,Ze),  geodetic  (ipg,Ag,hg),  and  ENU 
coordinate  systems  [7]. 


north  pole  (0°  latitude  90°  longitude),  and  the  Ye- axis  is  normal  to  the  XeZe- plane  in 
accordance  with  the  right  hand  rule. 

The  geodetic  frame  is  an  angular  coordinate  system  which  uses  an  ellipsoidal 
approximation  of  the  Earth  geoid  to  define  a  3D  point  in  terms  of  latitude  longitude 
(Ag),  and  altitude  (hg).  Latitude  is  the  angle  measured  from  the  equatorial  plane  to  the 
point  normal  to  the  surface  of  the  ellipsoid  and  ranges  from  -90°  to  90°.  Longitude  is  the 
angle  measured  counterclockwise  from  the  prime  meridian  in  the  equatorial  plane  and 
ranges  from  -180°  to  180°.  Altitude  is  the  height  above  the  surface  of  the  ellipsoid  along 
the  longitude  vector.  The  reference  ellipsoid  model  used  in  this  thesis  is  the  WGS  84 
ellipsoid  with  the  key  parameters  listed  in  Table  2.1. 

The  local  ECEF  reference  frame  is  specified  by  the  origin  located  at  a  point  p  in  ECEF 
coordinates  and  is  a  translation  of  the  ECEF  coordinate  system.  The  X'e,  Y',  and  Z'  axes  of 
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Table  2.1:  WGS  84  key  parameters  [8]. 


Semi-major  axis 

6378137  m 

Semi-minor  axis 

6356752.3142  m 

Eccentricity 

0.0818191908426 

Inverse  flattening 

298.257223563 

the  local  ECEF  reference  frame  are  parallel  to  the  Xe,  Ye,  and  Ze  axes  of  the  ECEF  frame, 
respectively. 

The  local  ENU  reference  frame  is  used  as  the  satellite  reference  frame  and  is  specified 
by  an  origin  and  the  vectors  east,  north,  and  up.  The  location  of  the  origin  is  located  at 
point  p  in  ECEF  coordinates.  The  up  vector  (U- axis)  points  along  the  longitude  vector 
and  is  normal  to  the  ellipsoid  surface.  The  north  vector  OV-axis)  points  towards  the  north 
pole  (Xe-axis)  and  is  tangent  to  the  ellipsoid  surface.  The  east  vector  (E-axis)  is  normal  to 
the  U N- -plane  in  accordance  with  the  right  hand  rule.  The  AE-plane  of  the  ENU  frame  is 
at  the  altitude  hg  of  p  and  is  tangent  to  the  surface  of  the  ellipsoid  at  Ag  and  tpg. 

The  sensor  reference  frame  is  used  to  specify  the  position  of  an  antenna  array  within 
the  local  satellite  ENU  reference  frame  and  consists  of  a  right  handed  orthogonal  set  of  x, 
y,  and  z  axes  originating  from  point  p.  The  relationship  between  the  sensor  frame  and 
ENU  frame  is  dependent  on  the  orientation  of  the  satellite  and  is  detailed  in  Section  3.2.  A 
summary  of  the  reference  frames  used  in  this  research  is  shown  in  Table  2.2  and  the 
angles  are  defined  in  the  following  sections. 

2.2  Multiple  Signal  Classification  with  a  Uniform  Circular  Array 

The  first  step  of  the  AOA  geolocation  process  is  to  determine  the  AOA  of  an  impinging 
RF  signal.  The  MUSIC  algorithm  with  a  Uniform  Circular  Array  (UCA)  antenna  utilizes 
the  phase  delays  of  an  impinging  Electromagnetic  (EM)  wavefront  across  the  individual 
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Table  2.2:  Reference  frames. 


Reference  Frame 

Origin 

Axes 

Angles 

ECEF 

Center  of  mass  of  the  earth 

XeYeZe 

Local  ECEF 

ECEF  position  p 

a,  s 

Local  ENU 

ECEF  position  p 

ENU 

4*  enu  >  &enu 

Sensor 

ECEF  position  p 

xyz 

antenna  elements  to  determine  the  AOA.  The  AOA  and  position  of  the  receiver  platform 
are  used  to  generate  a  3D  Line  of  Bearing  (LOB)  in  the  direction  of  the  emitter. 

2.2.1  Signal  Model. 

Consider  a  general  complex  baseband  signal  s(t)  of  unit  magnitude  \s(t)\  =  1  where 
Re{5(t)}  and  Im].v(/)}  are  the  In-phase  and  Quadrature-phase  components,  respectively. 
The  wireless  RF  transmission  of  s(t)  is  [9] 


st(t)  =  Re  I  ^Jojs(t)ej0)ct  j  (2.1) 

where  07  is  the  transmitted  signal  power  in  Watts,  a>c  =  2 nfc  is  the  carrier  frequency  in 
rad/s,  and  fc  is  the  carrier  frequency  in  Hz. 

The  signal  received  by  an  antenna  after  the  transmitted  EM  wave  propagates  through 
the  atmosphere  and  free  space  is 


sr(t)  -  Re  |  yJaj(t)s(t)eMtej(OD^ 


(2.2) 


where  cr2{t)  is  the  received  power  and  coD(t)  is  the  Doppler  frequency  shift  due  to  the 
relative  velocity  between  the  transmitter  and  receiver. 

The  received  power  is  given  as  [10,  11] 


crlit)  = 


<r2G,Gr 


l l J n  l ' s  (d(t),  /if) 


(2.3) 
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where  Gt  is  the  gain  of  the  transmit  antenna,  G,  is  the  gain  of  the  receive  antenna,  La  is  the 
atmospheric  loss,  Lp  is  the  polarization  mismatch  loss,  Ls  ( d(t ),  Ac )  is  the  free  space 
propagation  loss,  d(t )  is  the  Euclidean  distance  (range)  from  the  transmitter  to  receiver, 
and  Ac  is  the  wavelength  of  the  transmitted  signal.  The  free  space  propagation  loss  is  [10] 


where 


L,ldMM  =  ^j 


Ar  = 


fc 

and  c  =  299,792,458  m/s  is  the  speed  of  light  in  a  vacuum  [10]. 
The  Doppler  frequency  shift  is 


...  ~d(t) 

0Jn(t)  =  - coc 

c 


(2.4) 


(2.5) 


(2.6) 


where  d{t )  is  the  time  derivative  of  the  distance  (range  rate)  between  the  transmitter  and 
receiver.  Equation  (2.6)  is  derived  from  [5,  11] 


<jjr(t)  =  1 


d(t) 


U)c  —  COc  +  COo(t) 


(2.7) 


where  a>r(t)  =  2 nfr(t)  is  the  frequency  of  the  received  signal  in  rad/s  and  fr(t)  is  the 
received  frequency  in  Hz. 

The  received  signal  after  complex  demodulation  is  [9] 

s<i(t)  =  sr(t)e~J0Jct  =  yJcr2r(t)s(t)eJ0JD<t)t 
and  after  sampling,  with  a  period  of  Ts,  becomes 

sAn\  =  yj crl[n\s[n\e]0)DWn 


(2.8) 


(2.9) 


where  n  =  t\t=nTs  is  the  discrete  sample  index. 

Consider  a  planar  EM  wavefront  impinging  on  a  M  element  UCA  antenna  [12-16]. 
The  antenna  elements  lie  in  the  xy  plane  of  the  sensor  reference  frame  and  are  uniformly 
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spaced  around  the  circumference  of  a  circle  of  radius  r.  As  the  plane  wave  propagates 
across  the  UCA,  it  will  arrive  at  each  element  at  a  different  time  as  shown  in  Figure  2.2. 
The  reception  time  delays  are  observed  as  phase  delays  in  the  received  signals.  The  AOA 
of  the  wavefront  is  given  in  terms  of  elevation  and  azimuth  angles,  6  and  <p  respectively.  6 
is  measured  down  from  the  z  axis  and  0  is  measured  counterclockwise  from  the  x  axis  as 
shown  in  Figure  2.3.  The  observed  phase  delay  is  a  function  of  the  array  geometry  (radius 
and  element  spacing),  the  wavelength  of  the  impinging  signal,  and  the  AOA.  If  the 
receiving  antenna  elements  are  assumed  to  be  identical  and  omnidirectional,  the  phase 
delay  referenced  to  the  center  of  the  UCA  at  the  m-th  element  is  given  as  [15] 

am(9,  <p)  =  exp  r  sin(0)  cos  (0  -  ym)j  (2.10) 

where  A,-  =  c/ fr  is  the  wavelength  of  the  received  signal  and  ym  =  2n(m  -  1)  /M  is  the 
counterclockwise  angle  of  the  m-th  element  to  the  x  axis.  In  order  to  avoid  spatial  aliasing, 
adjacent  antenna  elements  are  required  to  be  spaced  no  further  apart  than  Ar/2  [12]. 


Figure  2.2:  Uniform  circular  array  phase  delay.  A  planar  wavefront  will  arrive  at  each 
antenna  element  at  a  different  time.  The  time  difference  is  observed  as  a  phase  delay  in  the 
received  signal  [9]. 
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Figure  2.3:  Uniform  circular  array  geometry.  M  identical  antenna  elements  lie  in  the  xy 
plane  of  the  sensor  reference  frame  and  are  evenly  spaced  along  the  circumference  of 
a  circle  of  radius  r.  For  an  impinging  planar  wavefront,  the  elevation  6  and  azimuth  (p 
angles  of  arrival  are  measured  down  from  the  z  axis  and  counterclockwise  from  the  x  axis, 
respectively  [16]. 


If  it  is  assumed  that  the  impinging  signal  is  narrowband  such  that  the  propagation  delay 
across  the  antenna  array  is  much  smaller  than  the  reciprocal  of  the  bandwidth  of  the 
impinging  signal  [9,  13],  then  the  signal  from  the  m-th  antenna  element  after 
demodulation  and  sampling  is  [15] 

xm[n\  =  am(6,  <f)sd[n ]  +  wm[n\  (2.1 1) 

where  wm[n\  is  a  zero  mean  complex  Additive  White  Gaussian  Noise  (AWGN)  random 
process  with  covariance  cr2I.  It  is  commonly  assumed  [12,  13,  16-18]  that  the  noise 
wm[n\  is  spatially  and  temporally  independent  of  the  signal  sd[n].  The  power  of  wm[n\  is 
given  as  [10,  11] 

=  xTsysW  (2.12) 

where  k  =  1.3806504  x  10~23  J/K  is  the  Boltzmann  constant  [10,  11],  Tsys  is  the  noise 
temperature  of  the  receiver  system,  and  W  is  the  bandwidth  of  interest. 
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Using  compact  matrix  notation  for  the  M  antenna  array  signals  with  one  impinging 


signal  present  and  N  samples, 


X  =  a  (0, 0)  s  +  W 


(2.13) 


where 


*i[l]  *i[2]  •••  Xl[N] 
x2[l]  x2[2]  •••  x2[N] 


=  x[l]  x[2]  •••  x[N]  (2.14) 

L  MxN 


xM[  1]  XM[2\  •••  xM[N] 


are  the  M  received  signals  as  defined  in  (2.11)  and  x[n]  is  the  received  signal  vector  at 


sample  n. 


a  (0,  (p ) 


«i(0,  (p) 
a2(0,  (p) 


(2.15) 


aM(0,<P) 

L  JMxl 

is  the  steering  vector  comprised  of  the  antenna  element  phase  delays  defined  in  (2.10), 


s=  Si[  1]  Sl[2]  •••  Sl[N] 


(2.16) 


is  the  sampled  demodulated  signal  and. 


Wi[l]  w  i[2]  •••  wi[A^] 

w2[l]  w2[2]  •  •  •  w2[iV] 


(2.17) 


wM[  1]  wm[2]  •••  wm[N ] 


are  independent  noise  realizations. 

The  general  signal  model  for  M  antenna  elements,  K  impinging  signals,  and  N  samples 


is  given  as  [9,  19] 


X  =  AS  +  W 


(2.18) 
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where  each  column  of  A  is  the  steering  vector  of  the  k-th  impinging  signal, 

A=  a(#i,0i)  a(02,02)  •••  <Pk)  (2-19) 

L  IMxK 

and  S  is  the  collection  of  K  demodulated  signals 

si[l]  sd2]  •••  sdN] 
s2[l]  s2[2]  •••  s2[N] 

sK[  1]  sK[  2]  •••  ,^[A] 

2.2.2  MUSIC  Algorithm. 

The  MUSIC  algorithm  [19]  utilizes  the  eigenstructure  of  the  spatial  covariance  matrix 
of  the  signals  received  by  M  antenna  elements  to  determine  the  AOA  of  K  impinging 
signals.  The  M  x  M  spatial  covariance  matrix  is  given  as  [9] 

R„  =  E  [XXH]  =  E  [(AS  +  W)  (AS  +  W)H]  =  A E  [SSH ]  A H  +  E  [wWff] .  (2.21) 

Since  it  is  assumed  that  the  noise  is  an  uncorrelated  zero  mean  complex  AWGN  process, 
the  noise  covariance  matrix  is 

E  [WW"]  =  ctIImxm  (2.22) 

and  (2.21)  becomes 

Rvx  =  ARSAAff  +  cr\lMxM  (2.23) 

where 

RSi  =  E  [SSH]  (2.24) 

is  the  K  x  K  signal  covariance  matrix.  In  the  case  of  zero  mean  uncorrelated  signals, 

R.s,  =  diag  [crj,  o\,  ■  ■  ■  ,  cr|]  (2.25) 

where  cr ],  cr\, . . . ,  cr|  are  the  signal  powers  [15,  20]. 
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Since  RTV  is  a  Hermitian  matrix  by  design  (Rvv  =  R"),  it  has  a  unitary 
eigendecomposition  with  real  eigenvalues  and  orthonormal  eigenvectors  [21].  When  the 
number  of  signals  is  less  than  the  number  of  antenna  array  elements  ( K  <  M )  the 
eigendecomposition  of  the  spatial  covariance  matrix  can  be  arranged  such  that  [9] 

„  [  if  D,  0  lr  f 

R»  =  QDQH=  Qs  Q„  Qv  Qw  (2  26) 

L  0  all  JL 

where  Q  v  consists  of  the  K  eigenvectors  which  span  the  signal  subspace  with 
corresponding  eigenvalues  Ds  =  diag  [di  +  oi, Ai  +  cr;„  •  •  •  ,AK  +  cr2w j  and  Qw  consists  of 
the  M  -  K  eigenvectors  which  span  the  noise  subspace  with  corresponding  eigenvalues 
cr\,.  The  form  of  (2.26)  comes  about  when  the  eigenvalues  contained  in  D  are  arranged  in 
descending  order. 

When  K  <  M  the  matrix  ARSSAH  is  singular  with  rank  K  and  the  relation  [19] 

det  (aR.s,A")  =  det  (Ra,  -  <r2j)  =  0  (2.27) 

indicates  that  cr2,  is  an  eigenvalue  of  Rn  which  occurs  M  -  K  times.  An  eigenvector  qM  of 
Rvr  corresponding  with  a  cr;,  eigenvalue  is  shown  to  be  orthogonal  to  the  columns  of  A  [9] 

R«qVv  =  (ARvvA//  +  oil)  qv,  =  0  +  cr^,Iqw  =  cr^,qlv  (2.28) 

and  the  M  -  K  qH,  eigenvectors  span  the  noise  subspace  defined  by  the  columns  of  QH  . 
There  are  K  linearly  independent  eigenvectors  q  v  with  corresponding  eigenvalues  d  v  of  the 
matrix  ARSSAH  which  are  also  eigenvectors  of  R„  as  shown  by  [9] 

Rxx-q,  =  (AR  vvA//  +  all)  q.y  =  d,q,  +  o^Iq,  =  (d,  +  oty  q,  (2.29) 

where  dA  +  o i  is  the  corresponding  eigenvalue  of  R  „.  The  K  q  v  eigenvectors  span  the 
signal  subspace  defined  by  the  columns  of  Qs.  Sorting  the  M  eigenvalues  of  Rvr  contained 
in  D  in  descending  order  (dj  +  cr;,  >  •  •  •  >  AK  +  cr;v  >  cr;,  >  •  •  •  >  oi), 

,  Ds  0 

D  =  diag  di  +oi,---  ,AK  +  ai,oi,---  ,oi  =  '  ,  (2.30) 

1  0  oi  I 

L  JMxM 
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and  arranging  the  corresponding  K  eigenvectors  of  Q  v  and  M  -  K  eigenvectors  of  Qu  to 
partition  the  Q  matrix 


Q  = 


Qv  Q, 


(2.31) 


yields  the  eigendecomposition  shown  in  (2.26).  Since  the  eigenvectors  contained  in  Q  are 
orthonormal,  the  signal  subspace  is  orthogonal  to  the  noise  subspace  (Qf  Q„,  =  0). 

The  orthogonality  of  the  signal  and  noise  subspaces  provides  a  means  to  determine  the 
AOA  of  the  impinging  signals.  The  columns  of  A  are  the  signal  steering  vectors  which  lie 
in  the  signal  subspace.  Due  to  subspace  orthogonality,  aH  ( 9 ,  cp)  Qw  =  0  when  9  and  cp  are 
the  angles  of  an  impinging  signal.  Since  the  steering  vector  is  a  function  of  the  known 
antenna  array  geometry,  evaluation  of  the  MUSIC  spectrum  [9,  19] 


MUSIC 


0 M ) 


(2.32) 


a"  (9,  cp)  QwQl!&  (9,  cp) 

over  the  range  of  9  and  (p  will  yield  peaks  whenever  6  and  (p  correspond  to  the  AOA  of  an 
impinging  signal.  The  location  of  the  peaks  are  taken  as  the  AOA  estimates  of  6  and  cp.  In 
practice,  Rxr  is  estimated  from  the  received  signal  data  [9,  12,  13] 


71=1 


(2.33) 


and  the  estimate  of  the  noise  subspace  Qw  is  used  in  (2.32)  with  a  grid  search  over  the 
range  of  9  and  cp.  An  example  of  the  MUSIC  spectrum  is  shown  in  Figure  2.4  with  peaks 
corresponding  to  the  AOA  of  the  impinging  signals. 


2.2.3  MUSIC  Theoretical  Performance. 

Analytic  expressions  for  the  theoretical  performance  of  MUSIC  using  a  UCA  have 
been  developed  in  [12,  13]  from  a  first  order  Taylor  series  expansion  of  the  MUSIC 
spectrum  in  the  vicinity  of  the  true  AOA.  The  expressions  are  asymptotic  for  a  large 
number  of  N  samples  and  yield  the  error,  covariance,  and  variances  of  the  estimated 
angles.  The  following  derivation  is  adapted  from  [12,  13]  in  terms  of  the  elevation  9  and 
azimuth  cp  angles  of  the  impinging  signals. 
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Figure  2.4:  Example  MUSIC  spectrum  for  a  4  element  UCA  with  2  impinging  signals.  The 
peaks  at  (6i,<f>i)  =  (60°,  120°)  and  {On,  (pi)  =  (25°,  270°)  correspond  to  the  AOAs  of  the 
impinging  signals. 

The  estimated  angle  error  vector  of  the  k- th  impinging  signal  in  the  sensor  reference 
frame  is  defined  as 


0k  ~  0k 


e e<!>,k  - 


(2.34) 


< Pk  -  (pk 


and  is  evaluated  as 


(2.35) 


where 


Re{ajQwQ%} 
Re(a"QM,.Q"a0}  a^Q.Q% 


(2.36) 


and 


-Reja^Q^Q^a,} 

-Re{a«QwQ%) 


(2.37) 
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The  terms  ae  and  a(j)  are  the  partial  derivatives  of  a  with  respect  to  6  and  cp. 


da 

a,~d0 


a  ,b  — 


da 

dcp 


jj-r  cos  ( 6 )  cos  (cp  -  yi)  exp  f  j  j-r  sin  (0)  cos  (cp  -  yi)] 

ffr cos  (6)  cos  (0  -  yM)  exp  [j^r  sin  (0)  cos  (cp  -  yM)\ 
-jj-/  sin  (0) sin  (0  -  7i)  exp  [j  Jr  sin  (0)  cos  (cp  -  yi)] 


(2.38) 


(2.39) 


-jfr  sin  ( 6 )  sin  (cp  -  yM )  exp  [jgr  sin  (0)  cos  (cp  -  yM)\ 

The  error  vector  of  the  fc-th  impinging  signal  is  Gaussian  distributed  with  zero  mean  and 
covariance  matrix 


Cnl  I,  — 


^Ip 

a?QwQ^ 

Re{a"Q„,Q*afl} 

2K  det  (E) 

Re{a?Q„Q%} 

J 

(2.40) 


e=ek,(j>=<pk 


where 


P  =  [R-]it  +  4R;;(A"A)-,R-' 


k,k 


(2.41) 


In  the  case  of  a  single  impinging  signal, 


p  =  crs2  +  crier  2  {aH a) 


-l 


cr:2. 


(2.42) 


2.2.4  CRLB  for  2D  Angle  Estimation  of  a  Single  Source  Using  a  UCA. 

This  section  derives  the  stochastic  Cramer-Rao  Lower  Bound  (CRLB)  for  2D  angle  of 
arrival  estimation  of  a  single  source  using  a  UCA  and  follows  a  similar  derivation  found 
in  [22] .  The  CRLB  is  the  lower  bound  on  the  variance  of  any  unbiased  estimator  of  6  and  cp 

var  cov  ( 6 ,  cp) 
cov  (cp,  6 )  var 

and  is  equal  to  the  inverse  of  the  Fisher  Information  Matrix  (FIM) 


>  CRLB  (e,  $) 


(2.43) 


CRLB  (e,f)  =  F 1-1  (0,0). 


(2.44) 
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From  (2.23)  the  received  signal  covariance  matrix  for  a  single  impinging  source  is 


Rvv  -  o-~aaH  +  ct'Imxm 


(2.45) 


where  a  is  the  steering  vector  as  defined  in  (2.15),  cr2s  is  the  signal  power,  and  cr2w  is  the 
noise  power.  Since  the  received  signal  vector  X  is  assumed  to  be  a  complex  valued 
Gaussian  distrusted  random  process  with  zero  mean  and  covariance  matrix  R„,  the 
stochastic  FIM  for  the  parameters  0  and  cp  is  given  as  [22,  23] 


F(0,0) 


trace 


- 1  dRx 

tx  38 

»-l<9Rv. 


trace  !  ffRx' If  R« 


trace 


trace 


<9R 


n-1  <9R.r.v  p 
38  -xx  3<j>  XX 


-1 


<9R 


_ 11-1  "Rr  1»  I 

3 0  ^XX  dif,  ^XX 


(2.46) 


where  the  partial  derivatives  of  the  covariance  matrix  with  respect  to  6  and  (p  are 


Element  1,1  of  the  FIM  is 


dR.vx 

dO 

dR.vx 

dep 


=  cr]aeaH  +  o^aaf 
=  cr]a^aH  +  aj  aaf. 


(2.47) 

(2.48) 


F i,i  (6,  cp)  =  trace  | ^pR«  ^^Rf' 


and  is  evaluated  by  substituting  (2.47)  into  (2.46) 


Fu  ( 0 ,  cp)  =  ((retrace  ((a.a^R^  +  aafR ~l)  (aeaffR;]  +  aafR^1)} 


Performing  the  matrix  trace  operation  yields 


F u  ( 6 ,  <p )  =  (o"j)2  (a//R“'ao)"  +  2  (a"R;>)  (af  R^a, e)  +  (af R^a)‘ 


The  terms  of  (2.51)  are  evaluated  by  making  use  of  the  following: 

2  m 

aHae  =  j—r  cos  ( 0 )  ^  cos  (<p  -  ym)  =  0 


m=  1 


af  a  =  -a^a^  =  0 


(2.49) 

(2.50) 

(2.51) 

(2.52) 

(2.53) 
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(2.54) 


H  2  n 

a%  =  j—r 

/li* 


M 


sin  (6)  ^  sin  (0  -  ym)  =  0 


m=  1 


a?a  =  -a^a,,,  =  0 


(2.55) 


since  for  a  UCA,  the  phase  center  of  the  array  is  located  at  the  center  of  the  array  [22] 
such  that 

M 


Z  cos  ( (p  -  ym)  =  0 

m=  1 
M 


Z  sin  (0  -  y,„)  =  0. 

m=  1 


(2.56) 


The  inverse  of  the  signal  covariance  matrix  is  evaluated  using  the  Woodbury  matrix 
identity  [22] 

(A  +  BCD r1  =  A-1  -  A~xb(c~x  +  DA~l  B^1  DA~l  (2.57) 


such  that 


r:!  =  -4 


-aa 


H 


+  aHa 


Evaluating  the  first  and  last  terms  of  (2.51)  yields 


a//Rja0  =  ~2 


a''a»-^,  ‘  „  (a"a)(aN 


+  a^a 


»"r;>  =  4- 


a^a 


A-(a"a)(a"a) 


4-  +  a^a 


=  0 


since  a?a  =  aHae  =  0.  The  second  term  of  (2.51)  is  evaluated  as 


a^Rja  =  -4 


aHa 


l—  (a"a)(a"a) 


+a»a 


and  after  simplification  with  a^a  =  ||a||z  and  a  common  denominator  becomes 


a//Rl'a  = 


1 


a 


^4  +  ||a||2‘ 

( T Z  11  11 


Similarly,  the  third  term  of  (2.51)  is  evaluated  as 

1 


a?R  la a  =  4- 

0  XX  V  rJ- 


afa0 


4+a"a 


(af  a)  (affa0) 


=  4||a,||2. 


(2.58) 

(2.59) 

(2.60) 

(2.61) 

(2.62) 

(2.63) 
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Continuing  with  the  evaluation  of  element  1,1  of  the  FIM,  (2.51)  becomes 


Fi,i  (6,(p)  =  2^-j— - - |N|2 

4  +  Hall2 


(2.64) 


and  with  the  assumption  of  omnidirectional  antenna  element  such  that  ||a||"  =  M  [22], 


(2.64)  becomes 


Evaluating  the  term 


Fu  (6,(p)  =  24-t^ — INI2- 
"4  +M 


(2.65) 


INI2  =  (2nr / Ar)2cos2  (6)  cos2  (6  -  ym) 


(2.66) 


and  noting  that 


cos2  (6  -  ym)  =  ^ 


1  +  cos  (26  -  2 ym)  M 


(2.67) 


yields 


and  (2.65)  becomes 


Nil2  =  — (2nrIAr)2 cos2  (0) 


Fu  (0, 4>)  =  4-j^ — (2^r / Ar)2cos2  (6) . 
W^  +  M 

<J% 


(2.68) 


(2.69) 


In  a  similar  manner,  elements  1,2  and  2,1  of  the  FIM  are  evaluated  as 


F i,2  (6,  (p)  =  F2,i  (6, 4>)  =  2^^- — (2 nrlAr)2  V  cos  (6  -  ym)  sin  (6  -  ym)  =  0  (2.70) 


ZM  '  M 

cos  (6  -  ym)  sin  (6  -  ym)  =)  \  sin  (26  -  2 ym)  =  0.  (2.71) 

m- 1  Z— im=l  ^ 


Element  2,2  of  the  FIM  is  evaluated  as 


F2>2  (6,(f>)  =  4^- — (2nrlAr)2sm-  (6) 
W^  +  M 

crz 


(2.72) 


sin2  (6  -  ym)  = 


1  -  cos  (26  -  2 ym)  M 


(2.73) 
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Using  (2.69),  (2.70),  and  (2.72)  the  FIM  for  the  parameters  0  and  0  is 


r2  M 2 


F«9,0)  = 


-(2  nr/Arf 


u-^  +  M 


cos2 ( 0 )  0 

0  sin2  (0) 


(2.74) 


and  the  AOA  CRLB  using  N  samples  of  the  signal  X  is  the  inverse  of  the  FIM  [22] 


CRLB  (0,0)  =  £[ F(0,0)]"l 


(2.75) 


Evaluating  (2.75)  with  (2.74)  produces  the  stochastic  CRLB  on  the  variance  of  2D  angle 
estimates  of  a  single  source  using  a  UCA  in  rad2 


+M 

<T% 


CRLB  (0, 0)  =  - 

4 NM2(2nr/Ar f 
<rw 

Defining  ij  as  the  Signal  to  Noise  Ratio  (SNR), 


1  /cos2  (0) 


0 


0 


1  /sin2  (0) 


(2.76) 


cr; 


0  = 


cri 


(2.77) 


the  elevation  angle  CRLB  is 


var 


(»)* 


i)  1  +  M 


t]N M2{2nr / Ar)2 cos2  (0) 


(2.78) 


and  the  azimuth  angle  CRLB  is 


var  0  > 


rj  1  +  M 


r]NM2(2nr/ d,  )2sin2  (0)' 

In  general,  the  CRLB  for  2D  angle  estimates  is  a  function  of  the  SNR  of  the  received 
signal,  number  of  samples  N,  number  of  antenna  elements  M ,  radius  of  the  UCA  r, 
wavelength  of  the  received  signal  Ar,  and  the  elevation  angle  0  of  the  received  signal. 


(2.79) 


2.3  Angle  of  Arrival  Geolocation 

Geolocation  of  a  RL  Signal  of  Interest  (SOI)  emitter  can  be  accomplished  with  a 
two-step  process  using  the  AOA  of  an  impinging  signal  and  the  position  of  the  receiver 
platform  when  the  signal  is  received.  The  first  step  in  the  process  is  to  determine  the  AOA 
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of  the  received  RF  signal  and  generate  a  LOB  in  the  direction  of  the  received  signal  from 
the  position  of  the  receiver  platform  when  the  signal  was  received.  As  the  receiver  passes 
over  the  emitter,  the  AOA  and  position  will  vary,  and  the  process  is  repeated  to  generate 
multiple  LOBs.  The  second  step  uses  the  collection  of  LOB s  to  estimate  the  position  of 
the  emitter. 

2.3.1  Lines  of  Bearing. 

Once  the  AOA  of  the  received  SOI  has  been  determined,  the  estimated  angles  along 
with  the  attitude  and  position  of  the  receiver  platform  are  used  to  generate  a  3D  LOB  in 
the  direction  of  the  emitter.  A  LOB  is  a  pointing  vector  originating  at  the  position  of  the 
receiver  p,  in  the  direction  of  the  estimated  angles  0,  and  (f)t.  The  attitude  of  the  receiver  is 
used  to  define  the  LOBs  in  a  common  coordinate  system  and  the  process  is  detailed  in 
Section  3.2.  In  the  general  case  shown  in  Figure  2.5,  the  LOB  vector  r,  from  p,  to  a 
geolocation  point  g  in  ECEF  coordinates  is  [24] 

r;  =  g  -  P,  (2.80) 

with  Cartesian  components 

n,x  =  gx  -  Pi,x 

ri,y=gy~Pi,y  ■  (2‘81) 

n,z  =  gz  ~  PUz 

The  azimuth  a,  and  elevation  s;  angles  in  the  local  ECEF  frame  are  [24] 

a,  =  atan2  (rUy,  ri  x) 

)  , - —  v  (2.82) 

£,  =  atari 2 y +  rfy,  ri  zJ 

where  a,  is  measured  counterclockwise  from  the  X'-axis  and  s,  is  measured  down  from 
the  Z'-axis.  The  angles  and  Cartesian  components  are  related  through 

n,x  =  llr.-ll  sin  (Si)  cos  (ad 

rity  =  INI  sin  (e,)  sin  (or,-)  •  (2.83) 

rUz  =  ||r,||  cos  (Sj) 
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The  vector  r,  is  the  true  (error  free)  LOB  from  p,  to  g  with  the  associated  bearing  angles 
a,  (g)  and  Si  (g).  The  estimated  z-th  LOB  in  the  ECEF  coordinate  system  used  for  the 
second  step  of  the  geolocation  process  is  defined  as  a  vector  originating  at  point  p,  and 
pointing  in  the  direction  of  d,  and  . 


Figure  2.5:  LOB  vector  r  pointing  from  p  to  g  with  angles  a  and  s  in  the  local  ECEF  frame. 


2.3.2  Angle  of  Arrival  Geolocation  Methods. 

The  second  step  of  the  AOA  geolocation  process  is  to  determine  the  position  of  a  RF 
SOI  emitter  from  the  LOB  data  sets  generated  in  the  first  step.  If  there  is  a  single  emitter 
or  multiple  emitters  that  have  been  segregated,  the  LOBs  as  shown  in  Figure  2.6  will  point 
in  the  direction  of  a  common  position  which  is  assumed  to  be  the  location  of  the  emitter. 
Two  methods  which  utilize  LOB  data  sets  to  generate  geolocation  estimates  are  the  Least 
Squares  (LS)  intersection  [25,  26]  and  Non-Linear  Optimization  (NLO)  [24,  27]. 

2.3.2. 1  Least  Squares  Intersection. 

Due  to  various  errors  (e.g.  AOA  estimation  error,  position  knowledge  error,  and 
attitude  knowledge  error)  the  collection  of  LOBs  will  not  all  intersect  at  the  location  of  the 
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Figure  2.6:  Three  lines  of  bearing  and  a  single  emitter. 

emitter.  However,  there  is  a  point  which  minimizes  the  squared  distances  between  that 
point  and  the  LOBs.  The  closest  point  is  found  as  the  LS  solution  that  minimizes  the  sum 
of  the  Euclidean  distances  squared  from  a  point  to  a  set  of  L  lines  [24].  The  LS 
intersection  point  and  distances  from  four  LOBs  is  shown  in  Ligure  2.7. 

The  method  used  in  [26]  to  find  the  LS  intersection  point  uses  the  estimated  bearing 
angles  d,-  and  st  of  the  z'-th  LOB  to  create  the  Cartesian  pointing  vector  from  p, 

sin  (£,)  cos  (ad 

u,  =  sin  (e,)  sin  (d,)  •  (2.84) 

COS  (£;) 

The  unit  vectors  of  the  L  LOBs  are  used  to  create  the  following  matrices 

L 

B  =  ^  1-3x3  -  Ujuf  (2.85) 

i=  1 
L 

b  =  (X3x3  -  u,uf)p;  (2.86) 

i=  1 

which  form  a  set  of  linear  equations 

Bg  =  b  (2.87) 

where  g  is  the  geolocation  estimate  of  the  emitter.  The  LS  minimum  distance  solution  [21] 
to  (2.87)  is 

gLS  =  (BrB)  “Brb  (2.88) 
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where  (b7  Bj  Br  is  the  Moore-Penrose  pseudoinverse  of  B  and  gL5  is  the  LS  intersection 
geolocation  estimate  of  the  emitter. 


Figure  2.7:  The  least  squares  intersection  is  the  point  which  is  the  closest  approximation 
to  the  intersection  of  a  set  of  lines  of  bearing.  The  location  of  the  point  is  determined  by 
minimizing  the  sum  of  the  distances  squared  to  all  lines  in  the  set  [24,  27]. 


2.3. 2.2  Non-Linear  Optimization. 

The  NLO  AOA  geolocation  method  as  described  in  [24,  27]  is  an  iterative  non-linear 
weighted  least  squares  approach  which  has  been  shown  to  produce  more  accurate 
estimates  of  an  emitter’s  location  than  the  LS  intersection  method.  The  variance  of  the 
angular  estimates  from  MUSIC  are  used  as  as  the  weighting  factors  for  the  collection  of 
LOBs.  NLO  also  produces  the  spatial  covariance  matrix  of  the  geolocation  estimate. 

The  iterative  NLO  process  begins  with  the  r,  vector  from  position  p,  to  a  geolocation 
point  g  as  defined  in  (2.80)  with  the  bearing  angles  at  and  e,  as  defined  in  (2.82).  The 
angles  or,-  and  e,  depend  on  the  non-linear  arctangent  function.  As  shown  in  Figure  2.8,  the 
arctangent  function  is  approximately  linear  with  unit  slope  in  the  range  of  -n/A  to  n/A.  In 
order  to  keep  a,  and  s,  in  the  linear  range  of  the  arctangent  function,  the  angles  are 
converted  to  one  of  four  local  NLO  reference  frames  shown  in  Figure  2.9.  The  frame 
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conversion  process  begins  by  using  the  following  rules  to  determine  which  reference 
frame  to  use  for  the  angle  pair: 

Frame  1:  if  -|  <  e,  <  |  and  <  at  <  |  or  ^  <  at  <  ^ 

Frame  2:  if  <  e,  <  |  and  |  <  a,  <  ^  or  <  a,-  < 

Frame  3:  if  |e,|  >  |  and  <  or,-  <  |  or  ^ 

Frame  4:  if  |e,|  >  |  and  |  <  a,-  <  ^  or  <  a,  <  -| 

Once  the  appropriate  reference  frame  has  been  identified,  the  angles  aNLO,i  and  sNLOj,  and 
the  gradients  with  respect  to  g,  Vg aNLOj  and  VgsNLOj  are  calculated  using  Equations 
(2.90)-(2.93)  for  the  selected  frame.  The  frame  conversion  process  ensures  the  values  of 
(Xnloj  and  snloj  He  in  the  -n/4  to  n/4  linear  range  of  the  arctangent  function  and  the 
appropriate  gradients  of  each  angle  pair  are  used  to  construct  the  Jacobian  matrix.  The 
following  term  is  defined  for  simplified  notation: 

||iT  J2  =  rUx2  +  rUy2.  (2.89) 


Figure  2.8:  Linear  approximation  of  the  arctangent  function  between  -n/ 4  and  nj4  [27]. 
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Figure  2.9:  The  four  NLO  reference  frames  [27]. 


Frame  1:  snloj  is  measured  from  the  X'Y'- plane  and  aNLo,i  is  measured  from  the 
X'-axis. 


sNLO,i  (g)  =  atan2  (ru,  Jr[x  +  r2vJ 
@NLO,i  (g)  =  atan2  (rUy,  rix) 

Vg Snloj  - 


V  g  C?NLO,i 


ri,xi\z 

r;H2||rw|| 

n,y 


a 


■wll 


(2.90) 


The  arctangent  function  with  two  arguments  (atan2)  is  used  to  return  angles  in  the  range 


of  -n  to  n  in  order  to  preserve  the  quadrant  of  the  angle. 
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Frame  2:  eNLOj  is  measured  from  the  X'Y'- plane  and  aNLOj  is  measured  from  the 


y'-axis. 


SNLOJ  (g)  =  atan2  f ru z,  ^ jr 2x  +  tj  ) 
&NLO,i  (g)  =  atan2  (rUx,  rUyJ 
V gSNLOJ  - 


V  g(*NLO,i 


(2.91) 


Frame  3:  sNLoj  is  measured  from  the  Z'-axis  and  aNLOj  is  measured  from  the  X'-axis. 


SNLO.i  (g)  =  atan2 


*1  + 


aNLOj  (g)  =  atan2  (rUy,  rUx) 

V gSNLOJ  ~ 


VgaNLOJ  - 


(2.92) 


Frame  4:  sNLOj  is  measured  from  the  Z'-axis  and  aNLOj  is  measured  from  the  y'-axis. 


&NLOJ  (g)  =  atan2  [  ^rfx  +  rfy,  ru 
aNLoj  (g)  =  atan2  (ritX,  rUy J 
Vg snloj  - 
^ga^Loj  - 


(2.93) 


0 


The  Jacobian  matrix  J  (g)  consists  of  the  partial  derivatives  of  sNLOj  and  a NLOj  with 
respect  to  the  x,  y,  and  z  components  of  g  and  is  given  as 


J  (g)  = 


dSNLO.l 

dgx 

dSNLO,  1 

dgy 

dSNLO,  1 
dgx 

VgS/VLO.l 

3«jvl0,i 

dgx 

daNLo,i 

dgy 

daNLO.l 

dgx 

VgCtNLOJ 

dSNLO.L 

dgx 

dSNLO.L 

dgy 

dSNLO,L 

dgx 

V g£NLO,L 

dciNLO,L 

L  dgx 

daNLo,L 

dgy 

damxo.L 
dgx  J 

2Lx3 

Vg  aNLO,L 

(2.94) 
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The  set  of  L  estimated  angle  pairs  from  MUSIC  in  the  appropriate  reference  frames  are 


contained  in  the  vector 


£1 


£nlo,  i 

&NLOA 

£nlo,l 

&NLO.L 


(2.95) 


J2Lxl 


and  the  angles  as  a  function  of  p,  to  g  in  the  corresponding  reference  frames  are 


£nlo,  1  (g) 

&NLOA  (g) 


«(g)  = 


£nlo,l  (g) 
O' NLO, L  (g) 


2Lxl 


(2.96) 


The  iterative  NLO  process  continues  by  using  the  LS  geolocation  estimate  as  the  initial 
NLO  estimate  (g0  =  g LS)  in  the  Gauss-Newton  iteration 


g/m  =  gA  +  (jjMJ^JjKAfi* 


(2.97) 


where 

An*  =  ft-ft(g*)  (2.98) 

is  the  difference  between  the  angles  estimated  from  MUSIC  and  the  angles  from  the  fc-th 
geolocation  estimate,  is  the  covariance  of  the  MUSIC  estimated  angles  which  are  used 
as  the  weighting  factors,  and 

Ja  =  J  (§a)  (2.99) 

is  the  Jacobian  of  the  Uth  geolocation  estimate.  Assuming  that  the  angles  estimated  from 
MUSIC  are  independent  and  Gaussian  distributed  implies  that  the  covariance  matrix 
of  the  angles  contains  variances  of  the  estimated  angles  along  the  diagonal. 
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Equation  (2.97)  is  iterated  until  the  difference  between  estimates  has  sufficiently 
converged  (||g*+i  -  g*||  <  e )  and  the  final  iteration  is  taken  as  the  NLO  geolocation 
estimate  of  the  emitter  (g nlo)- 

The  spatial  covariance  of  the  k-th  geolocation  estimate  is  related  to  the  angular 
covariance  matrix  through 

Sfc  =  (JKJ*)'1  (2.100) 

which  can  be  used  to  visualize  the  confidence  surface  using  the  method  developed  in 
Section  2.6.  The  NLO  method  will  be  shown  to  be  the  minimization  of  the  maximum 
likelihood  cost  function  in  Section  3.3  and  (2.100)  will  be  shown  to  be  equivalent  to  the 
Cramer-Rao  lower  bound  developed  in  Section  3.4. 


2.4  Maximum  Likelihood  Estimation,  Gauss-Newton  Iterations,  and  the 
Cramer-Rao  Lower  Bound 

Consider  a  general  known  scalar  nonlinear  function  /u  (6)  dependent  on  P  scalar 
parameters 


G  = 


(2.101) 


It  is  desired  to  estimate  the  P  unknown  parameters  from  N  >  P  noisy  observations  of  /j  ( 6 ) 
such  that 


x  =  n  (6)  +  w 


(2.102) 


where 


and 


T 


x  = 


p(0)  = 


[  Xi  ■■■  XN  j  , 
Pi(0)  •••  Pn(0) 


T 


(2.103) 

(2.104) 


w  = 


W 1  •  •  •  WN 


T 


(2.105) 
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The  noise  samples  w  are  assumed  to  be  independent  real  zero  mean  Gaussian  random 
processes  distributed  as 

w~JV(0,Cw)  (2.106) 


with  covariance  matrix 

Cw  =  diag [m;  j,  •  •  •  ,  (2.107) 

Due  to  the  noise  process,  the  N  observations  are  distributed  as 


x~N(n(0),  Cw) 


(2.108) 


with  the  conditional  Probability  Density  Function  (PDF) 

1 


/(x|0) 


:  exp 


1  iv 

-  E 


crt 


(2.109) 


V(2 n)N  det  (Cw)  v  -  „=i 

The  Maximum  Likelihood  Estimate  (MLE)  of  the  parameters  0MLE  is  taken  as  the  value  of 
0  which  maximizes  the  likelihood  function  [28] 


0Mle  =  argmax  {/  (x|  0)} 
6 


(2.110) 


and  is  equivalent  to  minimizing  the  term  in  the  exponent  of  the  PDF 


a  ■  ( V  (Xn  ~  (^))2 

Omle  =  arg  mm  1  ^ - v - 


.  n=  1 


o~z 


(2.111) 


Writing  (2. 1 1 1)  in  matrix  notation  yields 


Omle  =  arg  min  j[x  -^(0)]/'Cw1  [x-/i(0)]j.  (2.112) 

o 

Equation  (2.112)  is  a  nonlinear  weighted  least  squares  minimization  problem  which  can 
be  solved  using  the  Gauss-Newton  method. 

The  iterative  Gauss-Newton  method  uses  a  first  order  Taylor  series  expansion  to 
linearize  the  function  ju  (0)  about  the  value  0*  [28].  The  resulting  linear  least  squares 
problem  is  solved  and  the  produced  estimate  6k+{  is  used  to  linearize  n  (0)  for  the  next 
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iteration.  The  method  is  said  to  converge  to  a  solution  when  the  difference  between 
successive  estimates  is  sufficiently  small.  The  nonlinear  function  ft  (0)  is  approximated  as 
a  linear  function  about  0k  with  the  expansion 


n(0)*n(dk)  +  j(ok)  [o-ok] 


(2.113) 


where  the  Jacobian  matrix  of  partial  derivatives  of  n  ( 0 )  with  respect  to  the  P  unknown 
parameters  is 


3(0)  = 


<9/n(0) 

<90! 

00l(0) 

d6P 

V0P1  (0) 

dnN(0) 

L  a©! 

0yUjv(0) 
ddp  J 

NxP 

V, qPn  ( 0 ) 

(2.114) 


JVxP 


and  the  Jacobian  at  the  current  estimate  0k  is 

Ja=J(^).  (2.115) 

Substituting  the  linear  approximation  of  (2.113)  into  (2.112)  yields 

ek+ 1  =  argmin  {[x  -  (ft  (0k)  +  J,  [o  -  d*])f  C;1  [x  -  (ft  (ok)  +  3k  [0  -  0*])]}  (2.1 16) 

which  can  be  rearranged  as  [28,  29] 


0k+  i  =  argmin 
e 


lT 


x~  P  (@k)  +  3k0k  -3k0 


>1-1 


x~  P  (@k)  +  3k0k  -3k0 


=y 


0k+1  =  arg min  j[y  -  3k0\rCj  [y  -  3k0\\ 


(2.117) 

(2.118) 


where  all  the  terms  of  y  are  known  quantities.  Equation  (2.1 18)  is  a  weighted  linear  least 
squares  minimization  problem  of  0  with  the  solution  [28,  29] 


4,i  =  (JX'jT'jTCi'y- 

Back  substituting  for  y  produces 

oM  =  (jjc;1  j,)"1  jjc;1  [x  -  ft  (ok)  +  3 A] 


(2.119) 


(2.120) 
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which  after  distribution  of  terms 


4+i  =(j[c;1j4)“1jJc;1j,4  +  (j[c;1j,)“1j[c;1  [x-//(4)]  (2.121) 

.  . 

=1 

simplifies  to  the  Gauss-Newton  iteration  [28,  29] 

4+i  =  4  +  (jJc^j^j^c;1  [x-n  (4)1 .  (2.122) 

The  Gauss-Newton  method  is  said  to  converge  when  the  norm  of  the  difference  between 
successive  iteration  estimates  is  below  a  specified  threshold 

114+1 -4||<e  (2.123) 

and  the  final  iteration  is  taken  as  the  least  squares  estimate  of  the  parameter  vector 
4s e  =  4+i-  When  the  observations  are  Gaussian  distributed,  the  least  squares  estimate  is 
also  the  MLE  ( Omle  -  Qlse)  [28]. 

An  example  of  the  Gauss-Newton  iterative  minimization  method  for  2  unknown 
parameters  is  shown  in  Figure  2.10,  where  the  function 

f(0)  =  e\  +  6 102  +  4  (2.124) 

has  a  global  minimum  at  (0\,02)  =  (0, 0).  The  starting  point  for  the  first  iteration  step  is 
(0\ ,  02)  =  (-4,  -4)  and  the  norm  of  the  difference  from  the  20th  iteration  to  the  global 
minimum  is  less  than  10  4. 


The  CRLB  is  defined  [28]  as  the  statistical  lower  bound  on  the  variance  of  any 
unbiased  estimator  of  0  such  that 


var 


cov 


(ft) 


cov 


var 


(4,4) 

(ft) 


>  CRLB  (0)  (2.125) 


(4,4)  •• 

and  is  equal  to  the  inverse  of  the  FIM 

CRLB  (0)  =  F_1  (0) . 


(2.126) 
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Figure  2.10:  Example  Gauss-Newton  iterations  for  a  2  parameter  minimization  problem. 


Each  element  of  the  FIM  is  defined  as 


F  u(0) 


/5X(x|6/)\r/5X(x|  0)\ 

l  W  )  {  ddj  I 


where  the  log-likelihood  function  is  the  natural  logarithm  of  the  PDF 


(2.127) 


£(x|0)  =  ln(/(x|0)). 


(2.128) 


For  Gaussian  distributed  observations,  the  FIM  of  (2.127)  after  differentiation  and  the 
expectation  operation  becomes  [28] 


N 


Fy  (»)  =  X 

n=  1 


1  dnn  (0)  djj.n  (6) 

a-ln  ddi  dOj 


and  can  be  expressed  using  matrix  notation  as 


F  (0)  =  Jr  (0)  C“]  J  (0) . 


(2.129) 


(2.130) 


32 


For  a  large  enough  number  of  samples  N,  the  estimates  of  the  MLE  are  shown  to  be 
asymptotically  Gaussian  distributed  with  variance  equal  to  the  CRLB  such  that  [28] 

6~N(6,F-l(6)).  (2.131) 

2.5  Single  Sinusoidal  Signal  Parameter  Estimation 

This  section  derives  the  MLE  and  CRLB  for  estimating  the  parameters  of  a  single 
sinusoidal  signal  from  noisy  observations.  The  derivation  follows  the  example  found  in 
[28].  It  is  desired  to  estimate  the  normalized  frequency,  amplitude,  and  phase  of  a  single 
complex  exponential  signal  corrupted  by  AWGN.  Consider  N  discrete  samples  of  noisy 
observations  of  the  signal 


x  [/?]  =  crs  exp  [/  (2 nfn  +  (pf\  +  w  [ft] 


(2.132) 


where  the  noise  is  complex  Gaussian  distributed  as 

w[n]~c(0,a£l)  (2.133) 


and  the  signal  parameters  to  be  estimated  are  the  normalized  frequency  /,  amplitude  <xy, 
and  phase  (p 

a  =  [f,crs,<p]T.  (2.134) 

The  normalized  frequency  ranges  from  0  to  1  and  is  defined  as  f  =  fr/fs  where  fr  is  the 
received  frequency  of  the  signal  and  fs  is  the  sampling  frequency,  both  in  Hz.  Defining  the 
non-linear  function 

H  [ft]  =  cr s  exp  [j  (2 nfn  +  (p)]  (2.135) 


the  PDF  of  the  observed  signal  x  conditioned  on  the  parameters  a  is  [28] 


n- 1 


Ly  1/1  \  1 


^  ^  U  [ft]  -jj  \n\\2 

(2.136) 


w  n=  1 
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The  normalized  frequency  MLE  is  shown  [28,  30,  31]  to  be 


Imle  =  arg  max 


1  N 

Yj  x  M  exP  (rpxfn) 

n=  1 


(2.137) 


which  is  the  location  of  the  maximum  value  of  the  Discrete  Fourier  Transform  (DFT)  of  x 


fMLE  =  arg  max 


/ 


-DFT  [x] 

N 


(2.138) 


Similarity,  [28,  30]  the  amplitude  MFE  is  the  value  of  the  DFT  evaluated  at  fMLE 

1 


O'  s.MLE 


X 


DFT  [x] 


(2.139) 


f-ftALE 


An  example  of  the  frequency  and  amplitude  MFEs  is  shown  in  Figure  2.11. 


Figure  2.11:  The  normalized  frequency  and  amplitude  MFEs  of  a  single  complex 
exponential  signal  from  noisy  observations  are  the  location  of  the  maximum  value  and 
the  maximum  value  of  the  DFT  of  the  observations. 


The  FIM  for  complex  AWGN  is  given  as  [28] 


F, ,  («)  =  —Re 


cr- 


y  dnH  [«]  dfl  [ /?  ] 

2-f  da,  da ,• 

n=\  J 


(2.140) 
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and  the  partial  derivatives  of  n  \n]  with  respect  to  the  normalized  frequency,  amplitude, 
and  phase  are 


d/u  [/?] 


=  jlnncr s  exp  [j  (2 nfn  +  4>)\ 


dn  [n] 


d/u  [n\ 


=  exP  [j  (Znfn  +  (/>)] 
jcr s  exp  [j  (2nf. n  +  (p)] . 


Evaluating  (2.140)  yields 


F  (or)  = 


cr2s(2n)2  Yj  n2  0  cr227r  X  n 


N  0 


cr2s2n  X  n  0  Ncr2s 


Utilizing  the  following  summation  identities  [28], 


N(N+1) 


(2.141) 

(2.142) 

(2.143) 


(2.144) 


(2.145) 


N(N+  l)(2N+  1) 


(2.146) 


produces  the  CRLB  on  unbiased  estimates  of  the  single  sinusoidal  signal  parameters  [28] 


CRLB  (or)  =  F  1  (or)  = 


_ 6 _  0  _ — _ 

ri(2n)2N(N2-l)  T}2nN(N—\) 

0  ^  0 

J  2N  U 

-3  n  2N+1 

t]2nN(N-l)  V  tiN(N-  1) 


The  normalized  frequency,  amplitude,  and  phase  CRLBs  are  therefore 


(2.147) 


var(/)  > - y~ - - 

V'  rj{2n)2N{N2-  1) 

(2.148) 

2 

/A,  \  . 

var  (crs)  >  — 

(2.149) 

2iV  +  l 

var [(b)  >  - . 

~  ?]N  (N  -  1) 

(2.150) 

The  variance  of  the  received  frequency  fr  in  Hz2  is 

6/2 

van/,.  > - t - . 

V  ’  t](2tt)2N(N2  -  1) 

(2.151) 
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2.6  Confidence  Surfaces 


Using  a  single  unbiased  Gaussian  distributed  estimate  6  of  a  k  x  1  vector  parameter  6 
with  k  x  k  covariance  matrix  C§, 

e-N(e,Cfj)  (2.152) 

it  is  possible  to  define  a  Udimensional  surface  centered  at  the  estimate  0  which  has  the 
probability  P  of  containing  the  value  of  the  6  parameter.  The  ^-dimensional  confidence 
surface  is  defined  by  the  quadratic  relation 

[e  -  efc- 1  (e  -  e)  =  c2  (2.153) 


where  c  is  a  constant  to  be  determined  that  defines  a  contour  of  constant  probability.  The 
probability  P  that  the  value  6  is  contained  within  the  confidence  surface  centered  at  the 
estimate  6  is  given  as  [28,  32] 

P  =  Pr|(6/-0)rCl1(6/-0)  <  c2}.  (2.154) 

Letting  u  =  (o  -  6}  Cr 1  (()  -  0)  ,  since  0  is  Gaussian  distributed  with  covariance  Cg  and  u 
is  a  quadratic  sum  of  Gaussian  random  variables,  u  is  a  Chi-squared  random  variable  with 
k  degrees  of  freedom  [28,  32] 

u~xl  (2-155) 


The  PDF  of  the  Chi-squared  distribution  is  given  as  [33] 

xl  («) 


e~u/2uk/2-l 


2kl2Y(k/2) 

where  the  complete  Gamma  function  is  defined  as  [33] 


(2.156) 


oo 

J.-r-v, 


(2.157) 


Using  the  Cumulative  Distribution  Function  (CDF)  of  the  Chi-squared  distribution  to 
evaluate  (2.154),  the  probability  P,  constant  c,  and  dimension  k  are  related  through  [32] 

r  (IT) 


P  =  Pr  \u  <  c 


r(f) 


(2.158) 
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where  the  lower  incomplete  Gamma  function  is  defined  as  [33] 


y(a,[3)  =  J'  e  tf  ldt. 

o 

When  k  =  2  dimensions,  u  is  exponentially  distributed  and  the  constant  c  can  be 
defined  in  terms  of  the  probability  P  using  [28,  32] 


(2.159) 


(2.160) 


Using  c  and  (2.153),  the  resulting  surface  defines  a  confidence  ellipse  centered  at  6  with 
probability  P  of  containing  the  6  value  as  shown  in  Figure  2.12. 

When  k  =  3  dimensions,  a  value  of  c  ~  2.7959  with  (2.158)  yields  P  =  0.95,  and  is  a 
95%  confidence  ellipsoid  surface  defined  by 


(6/  -  dfc r1  (0-e)  =  (2.7959)2. 


(2.161) 


Figure  2.12:  2D  confidence  ellipse  centered  at  6  with  a  95%  probability  of  containing  Q. 
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III.  Methodology 


This  chapter  outlines  the  analysis  scenario  in  Section  3.1,  defines  the  relationship 
between  reference  frames  in  Section  3.2,  derives  the  geolocation  maximum  likelihood 
estimators  in  Section  3.3,  derives  the  Cramer- Rao  lower  bound  on  geolocation  estimates 
in  Section  3.4,  and  details  the  geolocation  performance  analysis  in  Section  3.5. 

3.1  Scenario  Overview 

Consider  a  single  RF  SOI  emitter  located  at  the  Air  Force  Institute  of  Technology 
(Latitude:  39.782°  N,  Longitude:  84.083°  W,  Altitude:  0  m)  transmitting  a  1315  MHz 
narrowband  signal.  The  signal  is  received  by  a  single  6U  CubeSat  in  a  60°  inclined  450 
km  circular  LEO  over  the  SOL  The  assumed  CubeSat  geolocation  payload  consists  of  a  4 
element  UCA  tuned  to  1315  MHz  and  calibrated  to  perform  AOA  measurements.  Each 
antenna  array  element  is  connected  to  a  RF  receiver  for  filtering,  demodulation,  and 
digitizing  the  received  signals  from  the  antennas.  The  4  payload  receivers  are  assumed  to 
be  frequency  and  phase  coherent  with  sufficient  bandwidth  and  sampling  rate  to  process 
the  received  signals.  Estimation  of  signal  parameters  from  the  received  signals,  running 
geolocation  algorithms,  and  interfacing  with  other  CubeSat  subsystems  is  accomplished 
through  a  dedicated  payload  processor  with  sufficient  processing  capability.  The  CubeSat 
is  also  assumed  to  have  a  Guidance  Navigation  and  Control  (GNC)  subsystem  consisting 
of  an  onboard  Global  Positioning  System  (GPS)  receiver  to  provide  position  and  velocity 
information,  and  an  Attitude  Determination  and  Control  System  (ADCS)  to  provide 
attitude  information  and  maintain  a  Local  Vertical  Local  Horizontal  (LVLH)  orientation 
over  the  RF  SOL  The  assumed  payload  functional  diagram  is  shown  in  Figure  3.1. 

This  analysis  scenario  considers  only  a  single  RF  emitter  with  no  co-channel 
interference.  It  is  assumed  that  the  signals  collected  from  multiple  RF  emitters  have  been 
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Uniform  Circular 
Antenna  Array 


Figure  3.1:  Assumed  CubeSat  geolocation  payload  consisting  of  a  4  element  UCA, 
4  frequency  and  phase  coherent  RF  receivers,  and  a  payload  processor  interfaced  with  the 
CubeSat  GNC  subsystem. 


properly  segregated  such  that  the  estimated  signal  parameters  used  as  inputs  into  the 
geolocation  algorithms  correspond  to  a  single  emitter.  Methods  for  data  association  and 
segregation  of  multiple  emitters  with  co-channel  interference,  as  well  as  the 
implementation  of  specific  RF  receiver  and  payload  processor  design,  are  beyond  the 
scope  of  this  thesis. 

The  payload  UCA  shown  in  Figure  3.2  consists  of  4  antenna  elements  evenly  spaced 
around  a  circle  on  the  nadir  face  of  a  6U  CubeSat.  The  UCA  elements  lie  in  the  .ry-planc 
and  are  centered  at  the  origin  of  the  sensor  reference  frame.  Using  (2.5)  the  wavelength  of 
the  1315  MHz  SOI  is  0.228  m,  and  the  diameter  of  the  UCA  is  0.114  m  to  satisfy  the  Tr/2 
antenna  element  spacing  requirement  of  MUSIC.  The  counterclockwise  angles  of  the 
antenna  elements  from  the  x-axis  are:  j\  -  45°,  y2  =  135°,  73  =  225°,  and  y4  =  315°.  It  is 
assumed  that  the  payload  antenna  array  is  properly  calibrated  with  a  uniform 
hemispherical  gain  pattern  to  receive  the  SOI  planar  EM  wavefront.  Implementation 
specific  antenna  design  to  include  variable  gain  patterns,  mutual  coupling,  calibration,  and 
other  antenna  parameters,  are  beyond  the  scope  of  this  thesis. 
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Figure  3.2:  Assumed  geolocation  payload  4  element  UCA  located  on  the  nadir  face  of  a 
6U  CubeSat. 


The  analysis  scenario  geometry  is  shown  in  Figures  3.3  and  3.4  where  the  RF  SOI  is 
located  at  position  g.  The  gain  pattern  of  the  RF  SOI  is  modeled  as  a  uniform  cone 
centered  at  g  and  extending  10°  above  the  horizon  with  an  80°  cone  half  angle.  The  single 
6U  CubeSat  receiver  platform  in  a  60°  inclined,  450  km  altitude,  circular  orbit  receives  the 
RF  SOI  at  ECEF  positions  p;  and  velocities  v,  while  traveling  through  the  cone.  Different 
types  of  RF  emitters  are  simulated  by  varying  the  number  of  signal  collects,  the  timing 
between  signal  collects,  and  the  SNR  of  the  received  signals.  Emitter  types  and  other 
simulation  variables  are  described  in  Section  3.5. 
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Figure  3.3:  Profile  view  of  the  analysis  scenario  geometry.  The  RF  SOI  gain  pattern  is  an 
80°  half  angle  cone  centered  at  g.  The  single  moving  platform  receives  the  SOI  at  positions 
p,  and  velocities  v,-. 


Figure  3.4:  Top-down  view  of  the  analysis  scenario  geometry. 
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3.2  Reference  Frame  Transformations 


The  relationships  between  the  local  reference  frames  used  in  this  thesis  are  shown  in 
Figure  3.5.  The  reference  frame  transformation  process  uses  the  points  p  and  g  with  the 
attitude  of  the  CubeSat  to  determine  the  azimuth  <p{ g)  and  elevation  0(g)  angles  in  the 
sensor  reference  frame  as  a  function  of  geolocation  point  g.  The  expressions  for  (f>( g)  and 
0(g)  are  used  to  derive  the  MLE  and  CRLB  for  AOA  geolocation  in  Sections  3.3  and  3.4. 


z 


Local  ECEF  Local  ENU  Sensor 

Figure  3.5:  Relationships  between  the  local  reference  frames. 


3.2.1  ECEF  to  Local  ECEF  Reference  Frame  Transformation. 

Using  the  LOB  vector  r  =  g  -  p  in  ECEF  coordinates  between  the  receiver  position  p 
and  geolocation  point  g,  the  corresponding  azimuth  and  elevation  angles  in  the  local 
ECEF  frame  are 

a  (g)  =  atan2  (gy  -  py,  gx  -  px) 
s  (g)  =  atan2  (  yj(gx  -  pxf  +  (gy  -  pyf,  gz  -  p. 


(3.1) 
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where  a(g)  is  measured  down  from  the  local  Z'-axis  and  s  is  measured  counterclockwise 
from  the  local  Z'-axis.  The  local  ECEF  LOB  vector  Cartesian  coordinates  are 


X'e  =  sin  (e  (g))  cos  (a  (g)) 
Y'e  =  sin  (e  (g))  sin  (a  (g)) 
Z'e  =  COS  (S  (g)) 


which  can  be  expressed  in  terms  of  p  and  g  as 


X'  = 


Y'  - 

1  e 


_ V( gx-Px)2+(gy-Py Y _ 

(gz-Pz)  1^-p^y-PyY^  /fezgd!  +7 

\  (gz-Pz)2  \  (Sx~Px)- 

_ (gy-Py)  a/( gx-Px)2+(gy~Py Y _ 

(gx-Px)(8z-Pz)  A/PS+i 


t.gx-px)2+(gy-py)2  |  ^ 
y  ( gz-Pz )2 

3.2.2  Local  ECEF  to  Local  ENU  Reference  Frame  Transformation. 

The  local  ENU  coordinates  are  obtained  through  the  Direction  Cosine  Matrix  (DCM) 
operation  on  the  local  ECEF  coordinates  [34] 

E  X'e 

N  =  R ENU/ECEFif’g^gj  Y'e  (3d 


where  the  DCM  Rrnu/ecef  dgj  is  a  function  of  the  geodetic  latitude  ipg  and  geodetic 
longitude  Ag  of  point  p, 

-  sin  (dg)  cos  (dg)  0 

Renu/ecef  [fs,  =  -  sin  (</?g)  cos  (dg)  -  sin  (</?g )  sin  (dg)  cos(</?g)  •  (3d 

cos  cos  (dgJ  cos  sin  (dg)  sin 
The  local  ENU  coordinates  in  terms  of  the  local  ECEF  coordinates  are 
E  =  -Xf  sin  (dg)  +  Yf  cos  (dg) 

N  =  -Xf  sin  (<^gj  cos  (dg)  -  Y'e  sin  f  </>gJ  sin  (dgJ  +  Z'e  cos  (</>g)  (3d 

U  =  X'e  COS  (<£g)  COS  (dg)  +  Y'e  COS  (</?g)  sin  (dg)  +  Z'e  sin  (<£g)  . 
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The  local  ENU  azimuth  angle  4>enu  (g)  is  measured  clockwise  from  the  77-axis  and  local 
elevation  angle  0ENU  (g)  is  measured  up  from  the  7/E-planc 


4* enu  (g)  =  atan2  (E,  N ) 

Oenu  (g)  =  atan2  (u,  V£2  +  TV2) . 


(3.7) 


The  expressions  for  4>enu  (g)  and  Oenu  (g)  in  terms  of  p  and  g  are  obtained  by  substituting 
(3.3)  and  (3.6)  into  (3.7). 

3.2.3  Local  ENU  to  Sensor  Reference  Frame  Transformation. 

It  is  assumed  that  the  CubeSat  is  in  a  LVLH  orientation  over  the  RF  SOI  where  the 
.v-axis  of  the  sensor  frame  is  aligned  with  the  velocity  vector  of  the  CubeSat  and  the  z-axis 
points  towards  nadir  in  the  opposite  direction  of  the  17-axis.  The  vy-planc  of  the  sensor 
frame  is  coplanar  with  the  TVE-plane  of  the  local  ENU  reference  frame.  The  relationship 
between  the  local  ENU  and  sensor  reference  frames  is  shown  in  Figure  3.6  and  the 
expressions  of  the  azimuth  (/>  (g)  and  elevation  0  (g)  angles  in  the  sensor  reference  frame 


3.2.4  Sensor  Reference  Frame  to  ECEF  Coordinate  System  Transformation. 

The  MUSIC  algorithm  described  in  Section  2.2  estimates  the  azimuth  and  elevation 
angles  of  the  received  SOI  which  are  used  to  generate  a  LOB  from  p  in  the  direction  of  the 

/V  A 

estimated  angles  f  and  6.  In  order  geolocate  the  RF  SOI,  the  LOB  is  transformed  from  the 
sensor  reference  frame  to  the  ECEF  coordinate  system.  The  process  begins  by 
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u 


Figure  3.6:  Relationship  between  the  local  ENU  and  sensor  reference  frames  assuming  a 
LVLH  CubeSat  orientation. 


transforming  the  estimated  angles  in  the  senor  frame  to  the  local  ENU  reference  frame, 


(pENU  —  <p  +  (pNorth  +  <P  error 
@ENU  @  2  ®error 

where  the  attitude  knowledge  error  is  expressed  in  terms  of  (perror  and  6error.  The  local 
ENU  Cartesian  components  are 


E  =  cos  (6enu)  sin  (< Penu ) 

N  =  cos  (Oenu)  cos  Penu )  (3.10) 

U  =  sin  (Oenu)  ■ 


The  local  ECEF  Cartesian  components  are  related  to  the  local  ENU  Cartesian  components 
through 


% 

Y'e 

-  Recef/enu 

E 

N 

_  %  _ 

U 

(3.11) 
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where  R ecef/enu  {pg,  4,,)  =  [R enu/ecef  (< Pg ,  /l,)]  is  the  DCM  from  the  local  ENU  to 
local  ECEF  reference  frame.  The  local  ECEF  estimated  azimuth  and  elevation  angles  are 


a  =  atan2  ( Y'e,X'e ) 

£  =  alan2(^)2  +  (^)2,2;j.  <3'’ 

The  /-th  estimated  LOB  in  the  ECEF  coordinate  system  is  the  vector  r,  originating  from 
the  position  p,  and  pointing  the  cr,  and  s,  direction, 


r  i  =  p  ,•  + 


K 


Z' 

i.e 


(3.13) 


where  p,  =  p,  +  perror  is  the  estimated  position  with  position  knowledge  error. 


3.3  Geolocation  Maximum  Likelihood  Estimators 

This  section  develops  the  geolocation  MLEs  used  in  this  thesis.  From  the  received 
signals,  the  collection  of  AOA  and  Frequency  of  Arrival  (FOA)  measurements  are  used  to 
estimate  the  location  g  of  the  RF  SOI.  From  (2.7),  the  received  FOA  in  Hz  consists  of  the 
RF  carrier  frequency  and  Doppler  shifted  frequency 

fr(g)  =  fc+fd(g).  (3.14) 


The  Doppler  frequency 


fAg)  =  ~~d(g)  (3.15) 

c 

is  a  function  of  the  position  g  of  the  SOI  along  with  the  position  and  velocity  of  the 
receiver  where  d  (g)  is  the  range  rate  of  the  distance  between  p  and  g.  Letting 

r  iT 


V  = 


VX  Vy  Vz 


be  the  XeYeZe  components  of  the  receiver  velocity  in  the  ECEF 


coordinate  system,  the  range  rate  is  given  as  [35,  36] 


dig)  = 


vrr(g) 

I|r(g)ll 


V,  iPx  ~  gx)  +  Vy  (Py  ~  gy)  +  C  iPz  ~  gz) 
y](Px  ~  gxf  +  ( Py  ~  gyf  +  iPz  ~  gz)2 


(3.16) 
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If  fc  is  assumed  to  be  known,  the  observed  frequency  after  demodulation  consists  solely  of 
the  Doppler  frequency  shift  and  the  FOA  of  the  i-th  signal  collect  is 

fc  Ev  ( PxJ  ~  gx)  +  Vy,i  (py,i  -  ,?v)  +  vzJ  ( pz,t  -  gz) 
fdji g)  =  --  - ln  2  '  (3'17) 

C  [  y](Px,i  -  gx f  +  (py,i  -  Sv)  +  ( Pz,i  -  gzf 

Consider  a  single  moving  receiver  in  the  ECEF  coordinate  system  at  position  p  and 
velocity  v  receiving  signals  from  a  single  RF  SOI  located  at  g.  At  each  signal  collect,  the 
angles  and  frequency  of  arrival  are  estimated  from  the  received  signal.  Let  the  i-th  set  of 
estimated  parameters  be 

0ENU,i 

Q  =  4>ENU,i  (3.18) 

hi 

where  6Enuj  and  (f)ENu,i  are  the  estimated  elevation  and  azimuth  angles  in  the  ENU 
reference  frame  to  account  for  attitude  knowledge  error,  and  hj  is  the  estimated  received 
frequency  of  the  signal.  The  true  parameter  values  as  a  function  of  p„  v,  and  g  are 

&ENU,i  (g) 

fii(g)=  <fevw(g)  (3-19) 

hi  (g) 

where  6ENuj  (g)  and  (pENuj  (g)  are  defined  in  (3.7),  and  fdi  (g)  is  defined  in  (3.17).  The 
variance  of  the  estimated  parameters  with  the  additional  attitude  and  frequency  knowledge 
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If  each  set  of  estimated  parameters  are  Gaussian  distributed  such  that 


ai~N(ai(g),Cci,l) 


(3.21) 


then  the  PDF  of  a  collection  of  L  sets  of  parameters  conditioned  on  g  is 

L  3  1  r  1  j 

f(Clu---,ClL  |g)  =  f[(2^)-2[det(Cn>i)]^exp  --(Q, -  Q, (g)f  .(g)) 


(3.22) 


and  taking  the  natural  logarithm  of  the  PDF  produces  the  log-likelihood  function 


£(Clu---  ,ClL |g)  = 

-  y  In  (2  tt)  -  l-  In  f\  det  (C^)  -  \  j]  (6,  -  Qf  (g)f  C^.  (4  -  Q,  (g)).  (3.23) 

L  /=1  J  /=  1 

If  the  parameter  estimates  are  independent  with  the  covariance  of  (3.20),  then  the 
log-likelihood  function  becomes 

£(Clu--- ,ClL  |g)  = 

3L  1  ^ 

-  —  In  (2n)  -  -  In  FI  (err  .  +  <T^ttitude ,)  (err.  +  cr^ttitude  (cri  +  cr|equency ,) 

L  (=i 

1  ^  ( 0ENU,i  ~  @ENU,i  (g))  (<f>ENU,i  ~  <pENU,i  (g))  ^  (fd,i  ~  fd,i  (g)) 

2  *  ^  (T2  4-  rr2  /  J  fy-  4-  rt -2  /  J  (J-2  4-  ct2 

I  1=1  OJ  attitude,/  /=  1  attitude,/  /=  1  frequency,/ 


(3.24) 


As  described  in  Section  2.4,  the  MLE  of  the  parameter  g  is  the  value  of  g  which 


maximizes  the  log-likelihood  function.  Maximization  of  the  log-likelihood  function  is 
accomplished  by  minimizing  the  the  summation  terms.  Utilizing  angle  and  frequency 
estimates,  the  AOA/FOA  geolocation  MLE  is  [37,  38] 


g AOA/FOA  ~ 

f  4-,  (OeNUJ  ~  OENU.i  (g))  4-,  (< t>ENU,i  ~  (pENU.i  (g))  4^  ( fd,i  ~  fd,i  (g)) 

\L  a2  +(r2 -  +  Zj  n-2  +(r2 -  +  L  - 

8  I  i=l  q  i  u  attitude,/  /=!  u  attitude,/  /=!  f  j  w  frequency,/ 


(3.25) 
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Using  just  the  estimated  angles,  the  AOA  geolocation  MLE  is  [39] 


(  Eh  ( @ENU,i  -  QeNUJ  (g))  ■A 

gA0A  =  argmin  n— r— 2 - +  Z, 


-  (pENU.i  (g)) 


— (j\  q-2  /  ^  q-2  _|_  0-2 

i=l  Qi  attitude,/  /=1  attitude,/ 


(3.26) 


Using  just  the  estimated  frequency,  the  FOA  geolocation  MLE  is  [5] 


(fdj  ~  fd,i  (g)) 

gFOA  =  arg  min  ^  -5— ~2 - 

8  I  ;=1  f  [  u  frequency,! 


(3.27) 


The  non-linear  weighted  least- squares  minimization  of  (3.25)  for  the  AOA/FOA  MLE 
is  solved  using  the  iterative  Gauss-Newton  process.  The  AOA  and  FOA  solutions  are 
derived  by  omitting  the  frequency  and  angle  terms,  respectively.  Equation  (3.25)  is 
expressed  in  matrix  notation  as 


g AOA/FOA  =  arg  min  {[Q  -  12  (gjfc  E  [ft  -  £1  (g)] 


(3.28) 


where  the  collection  of  L  sets  of  estimated  parameters  is 


A  A  p  A  A  A 

=  9eNU,  1  (pENU,  1  Jd,  1  ’  ’  ’  &ENU.L  (pENU.L  fd.L 


(3.29) 


with  covariance 


Cq  diag  Cq j  •  •  •  L 


and  the  true  parameter  values  as  a  function  of  g  are 


(3.30) 


r  13 

O  (g)  =  0ENU,  1  (g)  (pENU.i  (g)  fdA  (g)  -  -  -  0ENU.L  (g)  (pENU.L  (g)  fd.L  (g)  J 


(3.31) 


The  gradients  of  the  angle  and  frequency  functions  with  respect  to  the  ECEF  components 
of  g  are 


Vgifevt/,,  (g)  - 
^  gPpENU.i  (g)  = 


Vg fdj  (g)  = 


ddENUAg)  dSENU,i(g )  ddENUj(  g) 

dgx  dgy  dgz 

dlpENUjig )  9<!>ENuAg)  dlpENllA  g) 

dgx  dgy  dgz 


dfdA  g)  dfdAg)  dfdjig ) 

dgx  dgy  dgz 


(3.32) 


(3.33) 


(3.34) 
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The  partial  derivatives  of  the  non-linear  angle  and  frequency  functions  can  be  evaluated 
numerically  by  defining  the  following  perturbation  terms: 


1  0  0 


Av=  0  1  0 

A-=  0  0  1 


The  partial  derivatives  of  the  elevation  angles  are  approximated  as: 


(3.35) 


9SenuA%) 

dgx 

dSENU,i(  g) 
9gy 

d0ENU,i(g) 

dgz 


8ENU,i  (g  +  A  J  -  0ENU,i  (g) 
9 ENU,i  (g  +  A.)  -  Senuj  (g) 
9eNU,i  (g  +  Aj)  -  0ENU,i  (g)  • 


The  partial  derivatives  of  the  azimuth  angles  are  approximated  as: 


^£OT/.'(g)  _  (pENUi  (g  +  Ax)  -  4>ENU,i  (g) 
«  (pENUJ  (g  +  Ay)  -  (pENUJ  (g) 

dlpENUji  g) 


(pENU.i  (g  +  Aj)  -  (pENU,i  (g)  • 


(3.36) 


(3.37) 


The  partial  derivatives  of  the  frequencies  are  approximated  as: 

^7^  ~  fd,i  (g  +  Av)  -  fdj  (g) 

^OgT  ~  fdj  (8  +  Az)  -  fd,i  (g)  • 

The  gradients  of  the  z-th  set  of  true  parameters  is 

Vg9ENU,r  (g) 

VBa-(g)=  VgCpENuAS) 


VgfcU  (g) 


(3.38) 


(3.39) 


and  the  Jacobian  of  the  collection  of  L  sets  is 


Vgfli  (g) 


J(g)  = 


(3.40) 


Vg ClL  (g) 
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The  Gauss-Newton  iterations  for  the  AOA/FOA  geolocation  MLE  becomes 

g*+i  =  gk  +  (Jr  (at)  cA‘ J  (&))“'  f  (&)  c:1  [a -si  (§*)]  (3.4i) 

and  is  iterated  until  successive  estimates  converge  to  a  specified  threshold 

(llgi+i  -  g/t||  <  e)  or  a  specified  number  of  iterations  has  been  reached.  The  least  squares 

intersection  point  of  the  collection  of  LOBs  is  used  as  the  starting  iteration  (gt  =  g/s). 

The  Gauss-Newton  process  for  the  AOA  geolocation  MLE  is  equivalent  to  the  NLO 
method  described  in  Section  2. 3. 2. 2  and  [24,  27].  An  example  of  the  Gauss-Newton 
iterations  for  the  AOA/FOA  geolocation  MLE  is  shown  in  Figure  3.7  where  after  16 
iterations,  the  difference  between  iterations  is  less  than  1  m  (||gi6  -  gisll  <  1  m)  and  the 
final  iteration  is  the  AOA/FOA  estimate  of  the  emitter  location  (gAOA/FOA  =  gi6)- 


iteration  is  the  LS  intersection  of  the  collection  of  LOBs  and  the  final  iteration  is  the 
AOA/FOA  estimate  of  the  emitter  location. 
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The  95%  confidence  surface  of  the  final  geolocation  estimate  is  the  ellipsoid  in  the 


ECEF  coordinate  system  centered  at  g  and  defined  by 


/ 


x 


y 


V 


z 


T 

X 

X 

\ 

(jr(g)C^‘J(g))_l 

y 

-g 

/ 

z 

(2.7959)2 


(3.42) 


where  x,  y,  and  z  are  the  ECEF  components.  An  example  of  the  confidence  surfaces  for 
the  AOA,  FOA,  and  AOA/FOA  geolocation  estimates  is  shown  in  Figure  3.8  where  the 
term  (jr  (g)  Ct1  J  (g))  *  is  the  approximated  CRLB  for  estimate  g. 


Figure  3.8:  Confidence  surfaces  for  single  estimates  of  the  emitter  location  from  the  AOA, 
FOA,  and  AOA/FOA  maximum  likelihood  estimators.  The  FS  intersection  is  used  as  the 
initial  iteration  for  the  estimators. 
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3.4  Cramer- Rao  Lower  Bound  on  Geolocation  Estimates 


cov  (§)  >  CRLB  (g)  =  (jT  (g)  Cn]  J  (g))“‘  (3.43) 

For  a  particular  geolocation  analysis  scenario  where  L  sets  of  angle  and  frequency 
measurements  collected  by  a  single  moving  receiver  at  positions  p,  and  velocities  v,  from  a 
RF  emitter  located  at  g,  the  covariance  of  unbiased  estimates  of  g  is  bounded  by  the  CRLB 

var  (gx)  cov  (gx,  gy)  cov  (gx,  gz) 
co y(gx,gy)  var  (gv)  cov(gy,gz)  >  CRLB  (g) .  (3.44) 

co  v(gx,gz)  co  v(gy,gz)  var  (gz) 

The  FIM  for  Gaussian  distributed  AOA  and  FOA  measurements  is 

Faoa/foa  (g)  =  E  {[vgX (Gi ,  •  •  •  ,  ClL\  g)]  [VgX^Gi,  •  •  •  ,  Q/,|  g)j|  (3.45) 
and  after  evaluation  of  the  gradient  and  expectation  operations  becomes 

L 

F aoa/foa  (g)  =  [vga-  (g)]TCa)  [vga-  (g)].  (3.46) 

i=  1 

Using  the  matrix  notation  defined  in  Section  3.3,  the  FIM  for  AOA/FOA  is  expressed  as 

F  aoa/foa  (g)  =  Jr(g)Ca1J(g),  (3.47) 
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An  example  of  the  CRLBs  on  geolocation  estimates  is  shown  in  Figure  3.9  where  the 
CRLBs  are  visualized  as  95%  confidence  surfaces  centered  at  the  true  emitter  position  and 
defined  by  the  FIM  covariance  matrix.  Given  a  sufficient  number  of  L  signal  collects,  the 
estimates  of  the  emitter  location  are  asymptotically  distributed  as 

g  ~  JV  (g,  F"1  (g)) .  (3.51) 


Figure  3.9:  Cramer-Rao  lower  bounds  on  geolocation  estimates  of  the  emitter  location 
for  the  AOA,  FOA,  and  AOA/FOA  geolocation  methods  visualized  as  95%  confidence 
surfaces. 


3.5  Geolocation  Algorithm  Performance  Analysis 

This  section  describes  the  simulation  environment  used  to  evaluate  the  performance  of 
the  geolocation  algorithms  developed  in  this  thesis.  Systems  Tool  Kit  (STK)  10  developed 
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by  Analytical  Graphics  Inc.  was  used  to  model  the  analysis  scenario  geometry  of  a  single 
CubeSat  in  LEO  and  a  single  RF  SOI.  MATLAB  developed  by  Math  Works  Inc.  was  used 
implement  the  geolocation  algorithms  using  simulated  estimates  of  the  signal  parameters 
based  off  the  orbital  position  and  velocity  data  generated  from  STK. 

3.5.1  Orbital  Position  and  Velocity  Simulation. 

The  analysis  scenario  was  modeled  in  STK  and  outputs  the  orbital  position  and 
velocity  data  in  1  second  intervals  to  simulate  the  CubeSat  onboard  GNC  subsystem.  The 
RF  SOI  was  modeled  as  shown  in  Figures  3.3  and  3.4  as  an  80°  half  angle  cone  centered 
at  39.782°  N  latitude,  84.083°  W  longitude,  and  0  m  altitude.  The  single  CubeSat  platform 
was  modeled  as  a  60°  inclined,  450  km  circular  orbit.  The  right  ascension  of  the  ascending 
node  orbital  parameter  was  varied  in  2.5°  increments  to  produce  11  passes  of  the  CubeSat 
over  the  SOI  as  shown  in  Figure  3.10  with  the  time  durations  shown  in  Table  3.1.  Each 
pass  produces  a  set  of  ECEF  position  and  velocity  data  in  1  second  intervals  when  the 
CubeSat  is  traveling  through  the  SOI  cone.  The  intervals  of  pass  7  are  shown  in  Table  3.2 
as  an  example. 


Table  3.1:  Time  duration  of  each  orbital  pass  of  the  CubeSat  over  the  RF  SOI. 


Pass: 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

Duration  (s): 

72 

184 

247 

294 

331 

363 

386 

405 

417 

425 

428 

3.5.2  Geolocation  Algorithm  Simulation. 

The  analysis  scenario  positions  and  velocities  for  a  given  pass  are  imported  to  the 
MATFAB  environment  to  simulate  the  parameter  measurements  and  to  assess  the 
performance  of  the  geolocation  algorithms.  A  single  RF  SOI  is  located  at  position  g  and 
the  collection  of  L  measurements  are  simulated  from  the  true  parameters  of  the  pass  from 
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1 


11 


Figure  3.10:  11  orbital  passes  of  a  single  CubeSat  over  the  RF  SOL 


Table  3.2:  ECEF  position  and  velocity  data  for  pass  7. 


Time  t  (s) 

Xe(m) 

Ye  (m) 

Ze  (m) 

vx  (m/s) 

vv  (m/s) 

v,  (m/s) 

1 

-915252 

-5557270 

3860379 

4729.851 

2699.525 

5007.538 

2 

-910526 

-5554566 

3865384 

4731.386 

2705.758 

5002.701 

386 

968956 

-4096095 

5376484 

4934.411 

4771.814 

2746.139 

the  analysis  scenario.  All  measurements  are  assumed  to  be  Gaussian  distributed  with  the 
value  of  the  means  being  the  true  parameter  values.  The  measurement  variances  are 
defined  by  the  estimator’s  respective  CRLBs  and  an  additional  error  variance.  The  CRLB 
variances  are  in  terms  of  the  number  of  samples  N  per  measurement,  the  SNR  r\  of  each 
measurement,  and  the  associated  system  parameters.  The  performance  of  the  single 
sinusoidal  parameter  estimators  for  angle,  frequency,  and  SNR  measurements  will  be 
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shown  in  Section  4.1  and  the  threshold  is  defined  where  the  estimator  performance 
achieves  the  respective  CRLB. 


The  i-th  estimate  of  the  CubeSat  position  and  velocity  are  simulated  as 


and 


p,-~  W(p(-,mJositionI3X3) 


v,  ~  N  (v/,  ^elocilyI3x3) 


(3.52) 


(3.53) 


where  o-2  iMU(in  and  <x2elocity  are  variances  of  the  position  knowledge  error  and  velocity 
knowledge  error,  respectively. 

The  i-th  estimate  of  the  elevation  angle  in  the  local  ENU  reference  frame  is  simulated 
as 


Oenuj  ~  N  ( 0ENU,i  (g) ,  cr2g  j  +  cr2ttitude) 
where  the  true  elevation  angle  is 

Oeni hi  (g)  =  atan2  |t//,  ^Er+Nf^j , 


(3.54) 


(3.55) 


the  elevation  angle  CRLB  is 


CTn  ;  = 


[m  (g)]  1  +  M 


m  (g)  NM2{2nr/ArJf  cos2  ( 0ENuj  (g)  +  Jr/2)  (n/ 1 80)2  ’ 
and  cr2titude  is  the  variance  of  the  attitude  knowledge  error. 


(3.56) 


The  i-th  estimate  of  the  azimuth  angle  in  the  local  ENU  reference  frame  is  simulated  as 


<f>ENU,i  ~  N  ( <pENU,i  (g)  >  +  ^attitude) 


(3.57) 


where  the  true  azimuth  angle  is 


(Penuj  (g)  =  atan2  (Eh  N.) . 


(3.58) 


and  the  azimuth  angle  CRLB  is 


[>7i  (g)]  ‘  +  M 


oZ 


Vi  (g)  NM2(2nr/ArJf  sin2  (eENUJ  (g)  +  7t/2)  (tt/  1 80)2 


(3.59) 
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An  example  of  L  LOBs  in  the  ECEF  coordinate  system  using  the  estimated  positions  and 
angles  is  shown  in  Figure  3.11. 


Figure  3.11:  Fines  of  bearing  using  the  simulated  positions  and  angles  from  an  analysis 
scenario  pass. 

The  i-th  estimate  of  the  Doppler  shifted  frequency  is  simulated  as 


fd,i  ~  N  (fd,i  (g) ,  (r)d  i  +  cr^requency) 


(3.60) 


where  the  true  Doppler  shifted  frequency  is 


fdj  (g)  =  - 


(3.61) 


the  frequency  CRFB  is 


(3.62) 


rn(g)(2n)2N(N2-iy 


frequency 


is  the  variance  of  the  frequency  knowledge  error. 


The  /-th  estimate  of  the  received  SNR  is  simulated  as 


(3.63) 
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where  the  true  SNR  is 


m  (g)  = 


(jjGfGr 


LaLpLs. ,i  (g’  4)  KTSySW 


and  the  SNR  CRLB  is 


1 


cr„ 


(3.64) 


(3.65) 


v  2  N‘ 

The  received  SNR  is  a  function  of  the  assumed  transmitter  parameters,  assumed  receiver 
parameters,  and  free  space  path  loss.  If  the  transmitter  and  receiver  properties  are  assumed 
constant,  then  the  received  SNR  varies  due  to  the  free  space  path  loss,  such  that 

const. 


m  (g)  = 


where  the  free  space  path  loss  is  a  function  of  the  distance  between  transmitter  and 
receiver,  and  the  wavelength  of  the  transmitted  signal 


(3.66) 


f J.s./  (§,  dc)  — 


An 

Tr 


( PxJ  ~  gxf  +  (Py,i  ~  gyf  +  (PzJ  ~  gzf 


(3.67) 


In  order  to  generalize  the  received  SNR  parameter,  the  SNR  at  the  first  signal  collect  rj \  is 
used  to  define  the  SNR  of  subsequent  signal  collects  such  that 


n  i  (g)  =  m 


6 s, I  (§j  '1  r ) 

Ls,i  (g?  4-c) 


(3.68) 


An  example  link  budget  for  the  received  SNR  at  the  first  signal  collect  is  shown  in 
Table  3.3.  The  distance  from  transmitter  to  receiver  in  the  analysis  scenario  is 
approximately  1,600  km.  The  transmitter  parameters  are  expressed  as  the  Effective 
Isotropic  Radiated  Power  (EIRP).  The  atmospheric  attenuation  at  1,315  MHz  is  assumed 
to  be  0.5  dB  [40]  and  constant.  Assuming  a  vertically  polarized  transmit  antenna  and  a 
circularly  polarized  receive  antenna  array,  the  polarization  mismatch  loss  is  -3  dB. 
Assuming  a  noise  temperature  of  800  K  and  200  kHz  signal  bandwidth  of  interest,  the 
received  SNR  at  the  first  signal  collect  is  9.2  dB.  Since  defining  implementation  specific 
parameters  is  beyond  the  scope  of  this  thesis,  the  771  parameter  is  used  to  represent  the 
transmitter  and  receiver  system  parameters. 
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Table  3.3:  Example  link  budget  for  the  received  SNR  at  the  first  signal  collect. 


Transmitter  Frequency 

fc 

1,315 

MHz 

Propagation  Path  Fength 

d 

1,600 

km 

Transmitter  EIRP 

crjGt 

20.0 

dB 

Free  Space  Path  Foss 

Ls  (g,  Ar) 

-158.9 

dB 

Atmospheric  Foss 

La 

-0.5 

dB 

Receiver  Antenna  Gain 

Gr 

5.0 

dB 

Polarization  Foss 

LP 

-3.0 

dB 

Received  Signal  Power 

of 

-137.4 

dBW 

System  Noise  Temperature 

Tsys 

800 

K 

Noise  Bandwidth 

W 

200 

kHz 

Noise  Power 

vl 

-146.6 

dBW 

Received  SNR 

V 

9.2 

dB 

The  4  geolocation  algorithms  implemented  in  the  MATLAB  simulation  are  the  LS 
intersection,  and  the  AOA/FOA,  AOA,  and  FOA  MFEs.  Using  the  collection  of  L 
simulated  and  estimated  positions,  velocities,  angles,  frequencies,  and  variances,  the  FS 
intersection  algorithm  is  implemented  as  described  in  Section  2.3.2. 1  and  the  MFE 
algorithms  are  implemented  as  described  in  Section  3.3.  The  FS  intersection  estimate  is 
used  as  the  initial  estimate  (g!  =  gLS)  for  all  3  MFE  iterative  algorithms  and  the  estimated 
parameters  are  used  to  estimate  the  parameter  variances.  Examples  of  the  estimated  SNR, 
frequencies,  and  parameter  variances  for  39  signal  collects  are  compared  to  the  true  values 
in  Figures  3.12  through  3.16. 
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Doppler  Frequency  (kHz)  Received  SNR  (dB) 


Figure  3.12:  Estimated  received  SNR  per  signal  collect  for  771  =  0  dB. 


Figure  3.13:  Estimated  Doppler  frequency  per  signal  collect. 
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Azimuth  Variance  (deg )  Frequency  Variance  (Hz ) 


Figure  3.14:  Estimated  frequency  variance  per  signal  collect. 


x  10 


Figure  3.15:  Estimated  azimuth  variance  per  signal  collect. 


62 


Elevation  Variance  (deg) 


Figure  3.16:  Estimated  elevation  variance  per  signal  collect. 
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The  performance  of  the  4  geolocation  algorithms  is  assessed  by  comparing  the  Root 
Mean  Square  Error  (RMSE)  of  independent  Monte  Carlo  Simulation  (MCS)  trials  to  the 
RMSE  of  the  associated  CRLB.  An  example  of  200  MCS  geolocation  estimates  per 
algorithm  is  shown  in  Figure  3.17.  The  RMSE  of  the  CRLB  on  geolocation  estimates  is 
defined  as 

RMSEcslb  =  Vtrace{CRLB  (g)}.  (3.69) 


The  error  of  P  MCS  trials  is  defined  as 


Co  — 


llgi  -gll 


llgp  -gll 


JPxl 


The  sample  mean  of  the  MCS  trials  is 


/4 


-  1  V 

-  pZj  ezJ 


i=i 


and  the  sample  variance  is 

=  pTlE(*t  i-f')2- 

i=  1 

The  RMSE  of  the  MCS  trials  is 


(3.70) 


(3.71) 


(3.72) 


RMSEmcs  =  yjo-l+pl 


(3.73) 


3.5.3  Geolocation  Simulation  Parameters  and  Emitter  Types. 

This  section  defines  the  parameters  of  the  geolocation  simulation.  The  transmitter 
frequency  (fc  =  1,315  MHz),  number  of  antenna  elements  (M  =  4),  and  UCA  radius 
(r  =  AJ 4)  parameters  are  held  constant  throughout  the  simulations.  The  simulation 
parameters  listed  in  Table  3.4  are  varied  to  conduct  the  sensitivity  analysis  shown  in 
Section  4.2  to  assess  parameter  impact  on  the  accuracy  of  the  geolocation  algorithms.  The 
geometry  of  the  analysis  scenario  is  varied  through  the  1 1  orbital  passes  shown  in 
Figure  3.10  where  the  distance  from  transmitter  to  receiver,  relative  velocity,  angles, 
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Figure  3.17:  An  example  of  800  MCS  geolocation  estimates  (200  per  algorithm)  of  the 
position  of  a  RF  SOL  The  RMSE  of  the  MCS  estimates  is  compared  to  the  associated 
CRLB. 


frequencies,  and  total  pass  duration  are  unique  for  each  pass.  The  number  of  signal 
collects  parameter  L  is  varied  to  simulate  a  limited  number  of  signal  collects  during  a  pass 
where  the  collects  occur  sequentially.  For  example,  if  L  =  10  out  of  40  possible  signal 
collects,  the  signal  is  received  during  the  first  25%  of  the  pass.  The  number  of  samples  per 
collect  parameter  N  is  varied  to  simulate  different  signal  durations  and  sampling  rates. 

The  SNR  at  the  first  signal  collect  parameter  is  varied  to  simulate  different  transmitted 
signal  and  receiver  characteristics.  The  attitude  knowledge  error  parameter  cr^ttitude  is 
varied  to  simulate  the  error  from  the  CubeSat  ADCS  subsystem.  The  position  knowledge 
error  cr^osition  and  velocity  knowledge  error  cr“e]ot|ty  parameters  are  varied  to  simulate  the 
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error  from  the  CubeSat  GPS  subsystem.  The  frequency  knowledge  error  parameter 


^frequency *s  varied  to  simulate  local  frequency  oscillator  drift  and  carrier  frequency  offset 
errors. 


Table  3.4:  Geolocation  simulation  parameters. 


Orbital  pass 

p,  V 

Number  of  signal  collects 

L 

Number  of  samples  per  collect 

N 

SNR  at  the  first  signal  collect 

m 

Attitude  knowledge  error 

cr 2 
attitude 

Position  knowledge  error 

( T 2  •  ■ 
position 

Velocity  knowledge  error 

^"velocity 

Frequency  knowledge  error 

^"frequency 

Three  types  of  RF  SOI  emitters  are  simulated  by  varying  the  total  number  of  collects 
per  pass,  the  timing  between  collects,  the  number  of  samples  per  collect  and  the  received 
SNR.  A  high  power  spinning  radar  is  simulated  as  signal  collects  occurring  every  10 
seconds  with  a  relatively  high  SNR  and  a  relatively  low  number  of  samples  per  collect.  A 
moderate  power  burst  communications  emitter  is  simulated  as  signal  collects  occurring 
every  30  seconds  with  a  relatively  moderate  SNR  and  a  relatively  moderate  number  of 
samples  per  collect.  A  low  power  continuous  communications  emitter  is  simulated  as 
signal  collects  occurring  every  second  during  the  pass  with  a  relatively  low  SNR  and  a 
relatively  high  number  of  samples  per  collect.  The  3  emitter  types  are  summarized  in 
Table  3.5  and  the  performance  of  the  geolocation  algorithms  as  a  function  of  emitter  type 
is  compared  in  Section  4.2. 


66 


Table  3.5:  Simulation  parameters  for  three  emitter  types. 


Emitter  type 

Time  between 

signal  collects 

Received  SNR 

Number  of  samples 

Spinning  radar 

10  s 

High 

Low 

Burst  communications 

30  s 

Moderate 

Moderate 

Continuous  communications 

1  s 

Low 

High 
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IV.  Results  and  Analysis 


This  chapter  reports  the  statistical  performance  of  the  single  sinusoidal  signal 
parameter  estimators  and  the  geolocation  algorithms  used  in  this  thesis.  Section  4. 1 
reports  the  performance  of  the  angle,  frequency,  and  amplitude  estimators  in  terms  of  the 
SNR  and  number  of  samples,  and  defines  the  performance  threshold  where  the  estimator 
variance  achieves  the  respective  CRLB.  Section  4.2  reports  the  impact  of  varying  the 
system  parameters  on  the  geolocation  accuracy  of  the  geolocation  algorithms  through  a 
parameter  sensitivity  analysis. 

4.1  Single  Sinusoidal  Signal  Parameter  Estimator  Performance 

This  section  evaluates  the  performance  of  the  angle  (MUSIC),  frequency  (MLE),  and 
amplitude  (MLE)  parameter  estimators  used  in  this  thesis  as  compared  to  the  respective 
CRLBs.  The  threshold  where  the  estimator  variances  achieve  the  CRLB  was  determined 
in  terms  of  the  SNR  r/  and  number  of  samples  N  of  the  received  signal.  Since  it  is  assumed 
that  there  is  only  a  single  emitter  present  with  no  co-channel  interference,  the  received 
signal  is  modeled  as  a  single  complex  exponential  signal  consisting  of  N  samples  with 
normalized  frequency  /  at  a  SNR  of  tj  =  cr2s  such  that 

5  \n]  =  crs  exp  [jlnfn] .  (4.1) 


The  signal  is  phase  shifted  according  to  (2.10)  with  a  simulated  4  element  UCA  with 
radius  r  =  Ar/ 4  and  true  angles  of  9  =  n/ 4  and  cf>  -  7 r/4 


exp(jf  sin  (0)  cos  i 

*-ij) 

exp(j'|  sin  (0)cos(<; 

»-¥)) 

exp(j'|  sin  (0)cos(<; 

»-?)) 

exp(j'|  sin  (0)cos(<; 

(4.2) 
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Each  of  the  4  phase  shifted  received  signals  is  corrupted  with  independent  unit  power 
(cr2  =  1)  complex  AWGN  such  that 

xm[n\  =  am(6,(p)s[n ]  +  wm[n ]  (4.3) 

and  the  4  x  N  samples  of  the  simulated  received  signals  are  contained  in 

X  =  a(0, 0)s  +  W.  (4.4) 


The  AOA  of  the  received  signals  are  estimated  from  the  estimated  spatial  covariance 
matrix  R„  =  -f  XX;/  and  a  2D  angle  grid  search  over  the  MUSIC  spectrum 


1 


=  arg  max  • 

9,<t> 


aH  ( 6 ,  cp)  Q„  Q^a  (0,  <p) 


(4.5) 


The  location  of  the  peaks  of  the  MUSIC  spectrum  are  taken  as  the  AOA  estimates  of  the 
received  signal.  The  normalized  frequency  and  amplitude  are  estimated  from  the  DFT  of 
the  first  received  signal  and  the  SNR  is  estimated  as  r)  -  cr2 


f  =  arg  max 

/ 


^DFTM 


<x2  = 


^DFT  [xj 


2 

/=/ 


(4.6) 

(4.7) 


The  associated  angle,  frequency,  and  SNR  CRFBs  as  a  function  of  the  SNR  and  number 
of  samples  of  the  received  signals  are 


CRFB 


CRFB  (/) 


?/  1  +  4 


/7A42(;r/2)2sin2  (n/4)  (tt/180)2 
6 


CRLB  (fi)  = 


rj  (A3  -  N)  {In)2 

1 


(4.8) 

(4.9) 
(4.10) 


The  following  figures  report  the  variance  of  the  angle,  frequency,  and  SNR  estimators  of 
noisy  signals  over  a  range  of  N  and  //  values.  Each  variance  is  calculated  from  3,000 
independent  MCS  trials  per  N  and  rj  value,  where  N  ranges  from  10  to  1,000  samples  and 
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tj  ranges  from  -20  to  30  dB.  At  each  grid  point,  the  associated  CRLB  is  calculated  for 
performance  comparison. 

The  surface  plot  shown  in  Figure  4. 1  reports  the  performance  of  the  frequency  MLE 
compared  to  the  frequency  CRLB  over  the  range  of  N  and  rj  values.  At  low  values  of  N 
and  T],  the  MLE  variance  approaches  a  value  ofl2/12  =  8.33x  10-2,  which  is  the  variance 
of  a  uniform  random  variable  ranging  from  0  to  1.  The  values  of  N  and  //  where  the 
variance  of  the  frequency  MLE  approaches  the  CRLB  is  defined  as  the  frequency 
estimator  performance  threshold.  Cross  sections  of  Figure  4. 1  at  constant  values  of 
N  =  170  samples  and  rj  =  —  4  dB  are  shown  in  Figures  4.2  and  4.3  respectively. 


Figure  4.1:  Performance  of  the  frequency  MLE  compared  to  the  frequency  CRLB  as  a 
function  of  SNR  and  N. 
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Figure  4.2:  Performance  of  the  frequency  MLE  compared  to  the  frequency  CRLB  as  a 
function  of  SNR  with  constant  N  =  170  samples. 


Figure  4.3:  Performance  of  the  frequency  MLE  compared  to  the  frequency  CRLB  as  a 
function  of  N  with  constant  SNR  =  -4  dB. 
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The  surface  plot  shown  in  Figure  4.4  reports  the  performance  of  the  azimuth  angle 
MUSIC  estimator  compared  to  the  azimuth  angle  CRLB  over  the  range  of  N  and  //  values. 
At  low  values  of  N  and  rj,  the  MUSIC  variance  approaches  a  value  of 
3602/ 12  =  1.08  x  104,  which  is  the  variance  of  a  uniform  random  variable  ranging  from  0 
to  360.  The  MUSIC  variance  is  less  than  the  angle  CRLB  at  low  values  of  N  and  r\ 
because  the  angle  CRLB  assumes  a  Gaussian  distribution  rather  than  the  uniform 
distribution.  The  values  of  N  and  rj  where  the  variance  of  the  MUSIC  angle  estimator 
approaches  the  CRLB  is  defined  as  the  angle  estimator  performance  threshold.  Cross 
sections  of  Figure  4.4  at  constant  values  of  N  =  20  samples  and  rj  -  -10  dB  are  shown  in 
Figures  4.5  and  4.6  respectively. 


SNR  (dB) 


Figure  4.4:  Performance  of  the  MUSIC  azimuth  angle  estimator  compared  to  the  azimuth 
angle  CRLB  as  a  function  of  SNR  and  N. 
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Figure  4.5:  Performance  of  the  MUSIC  azimuth  angle  estimator  compared  to  the  azimuth 
angle  CRLB  as  a  function  of  SNR  with  constant  N  =  20  samples. 


Figure  4.6:  Performance  of  the  MUSIC  azimuth  angle  estimator  compared  to  the  azimuth 
angle  CRLB  as  a  function  of  N  with  constant  SNR  =  -10  dB. 
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The  surface  plot  shown  in  Figure  4.7  reports  the  performance  of  the  SNR  MLE 
compared  to  the  SNR  CRLB  over  the  range  of  N  and  t]  values.  The  SNR  MLE  variance  is 
less  than  the  SNR  CRLB  at  low  values  of  N  and  //  because  the  SNR  CRLB  is  only  valid 
for  unbiased  estimates.  The  values  of  N  and  rj  where  the  variance  of  the  SNR  MLE 
approaches  the  CRLB  is  defined  as  the  SNR  estimator  performance  threshold.  Cross 
sections  of  Ligure  4.7  at  constant  values  of  iV  =  170  samples  and  rj  =  -10  dB  are  shown  in 
Ligures  4.8  and  4.9  respectively. 
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Ligure  4.7:  Performance  of  the  SNR  MLE  compared  to  the  SNR  CRLB  as  a  function  of 
SNR  and  N. 
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Figure  4.8:  Performance  of  the  SNR  MLE  compared  to  the  SNR  CRLB  as  a  function  of 
SNR  with  constant  N  =  170  samples. 


Figure  4.9:  Performance  of  the  SNR  MLE  compared  to  the  SNR  CRLB  as  a  function  of  N 
with  constant  SNR  =  -10  dB. 
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The  ratios  of  the  estimator  variances  to  the  corresponding  CRLBs  of  Figures  4.1,  4.4, 
and  4.7  are  expressed  in  dB  (101ogK)  [variance/CRLB])  and  are  shown  in  Figures  4.10 
through  4.12  for  the  frequency,  angle,  and  SNR  estimators  performance,  respectively.  The 
+  1  dB  contours  of  the  ratio  plots  are  shown  in  Figure  4.13  where  the  value  of  the 
estimator  variance  is  within  1  dB  of  the  associated  CRLB.  The  frequency  contour  plot  is 
similar  to  [41].  As  long  as  the  values  of  N  and  rj  are  within  the  region  to  the  right  of  the 
frequency  1  dB  threshold  contour,  the  use  of  Gaussian  distributed  signal  parameter 
estimates  in  the  geolocation  algorithm  simulation  with  the  variances  defined  by  the 
associated  CRLB  as  described  in  Section  3.5.2  is  justified. 


Figure  4.10:  Ratio  of  the  frequency  MLE  variance  to  the  frequency  CRLB,  expressed  in 
dB,  over  the  range  of  N  and  SNR. 
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Figure  4.11:  Ratio  of  the  MUSIC  azimuth  angle  estimator  variance  to  the  angle  CRLB, 
expressed  in  dB,  over  the  range  of  N  and  SNR. 


Figure  4.12:  Ratio  of  the  SNR  MLE  variance  to  the  SNR  CRLB,  expressed  in  dB,  over  the 
range  of  N  and  SNR. 
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Figure  4.13:  Contour  plot  of  the  frequency,  angle,  and  SNR  estimator  performance  ratios 
showing  the  +1  dB  thresholds. 
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4.2  Geolocation  Algorithm  Performance  Sensitivity  Analysis 

This  section  reports  the  impact  of  various  system  parameter  values  on  the  RMSE  of  the 
4  geolocation  algorithms  used  in  this  thesis.  The  system  parameters  include:  the  orbital 
pass  used  in  the  analysis  scenario,  number  of  signal  collects  during  the  pass,  number  of 
samples  per  signal  collect,  SNR  at  the  first  signal  collect,  frequency  knowledge  error, 
attitude  knowledge  error,  position  knowledge  error,  and  velocity  knowledge  error.  The 
spinning  radar,  burst  communications,  and  continuous  communications  emitter  types  are 
used  for  the  sensitivity  analysis  which  is  implemented  as  the  geolocation  simulation 
described  in  Section  3.5.2.  The  following  default  system  parameters  are  used  for  all  3 
emitter  types:  pass  number  7,  5  Hz  frequency  knowledge  error,  0.2  deg  attitude 
knowledge  error,  10  m  position  knowledge  error,  and  1  m/s  velocity  knowledge  error.  A 
sampling  frequency  of  400  kHz  (Ts  =  2.5  [ is)  is  used  for  the  simulated  frequency 
estimates.  The  system  parameter  under  consideration  is  varied  while  the  others  are  held 
fixed  for  each  scenario.  The  RMSE  of  the  4  geolocation  algorithms  is  calculated  from 
3,000  independent  MCS  trials  per  point  and  compared  to  the  corresponding  CRLB. 

The  spinning  radar  emitter  considered  in  Section  4.2.1  is  implemented  with  the 
following  parameters:  20  dB  SNR  at  the  first  signal  collect,  100  samples  per  signal 
collect,  and  10  seconds  between  signal  collects. 

The  burst  communications  emitter  considered  in  Section  4.2.2  is  implemented  with  the 
following  parameters:  10  dB  SNR  at  the  first  signal  collect,  400  samples  per  signal 
collect,  and  30  seconds  between  signal  collects. 

The  continuous  communications  emitter  considered  in  Section  4.2.3  is  implemented 
with  the  following  parameters:  0  dB  SNR  at  the  first  signal  collect,  800  samples  per  signal 
collect,  and  1  second  between  signal  collects. 

The  geolocation  performance  of  the  AOA/FOA  geolocation  MLE  for  the  3  emitter 
types  is  compared  in  Section  4.2.4. 
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4.2.1  Spinning  Radar  Emitter  Sensitivity  Analysis. 

Figure  4.14  reports  the  RMSE  of  the  4  geolocation  algorithms  and  associated  CRLBs 
for  the  1 1  orbital  passes  in  the  analysis  scenario.  The  orbital  pass  impacts  the  position, 
velocity,  and  total  number  of  signal  collects  used  by  the  geolocation  algorithms.  The  AOA 
LS  intersection  has  similar  performance  to  the  AOA  MLE  for  passes  1-4  due  to  similar 
valued  angle  measurement  variances  and  relatively  low  number  of  total  signal  collects.  As 
the  number  of  signal  collects  increases  with  a  wider  range  of  estimate  variances  for 
weighting  of  the  angle  measurements,  the  AOA  MLE  has  lower  RMSE  than  the  AOA  LS. 
The  AOA/FOA  MLE  geolocation  algorithm  has  a  lower  RMSE  than  the  other  3 
algorithms  for  all  passes. 


□ 

AOA  LS 

« 

AOA  MLE 

o 

FOA  MLE 

0 

AOA/FOA  MLE 

AOA  CRLB 

— 

FOA  CRLB 

AOA/FOA  CRLB 

Figure  4.14:  Sensitivity  analysis  of  varying  the  orbital  passes  with  the  spinning  radar 
emitter. 
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Figure  4.15  reports  the  RMSE  as  the  total  number  of  signal  collects  out  of  a  possible  39 
total  signal  collects  for  pass  number  7  is  varied.  The  RMSE  of  all  4  algorithms  decreases 
as  the  number  of  signal  collects  increases  and  begins  to  level  out  around  25/39  signal 
collects.  The  AOA/FOA  MLE  geolocation  algorithm  has  the  lowest  RMSE  throughout. 
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Figure  4.15:  Sensitivity  analysis  of  varying  the  number  of  collects  with  the  spinning  radar 
emitter. 
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Figure  4.16  reports  the  RMSE  as  the  number  of  samples  per  signal  collect  is  varied. 
The  RMSE  decreases  as  the  number  of  samples  per  signal  collect  increases  until  the 
frequency,  attitude,  position,  and  velocity  knowledge  errors  determine  the  performance  of 
the  algorithms,  at  which  point  the  RMSE  begins  to  level  out.  The  RMSE  of  the  FOA  MLE 
geolocation  algorithm  decreases  at  a  greater  rate  due  to  the  1  / N 3  term  in  frequency  CRLB. 
The  AOA/FOA  MLE  geolocation  algorithm  has  the  lowest  RMSE  throughout. 
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AOA  LS 
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AOA  MLE 
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FOA  MLE 
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AOA  CRLB 
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FOA  CRLB 

AOA/FOA  CRLB 

Figure  4.16:  Sensitivity  analysis  of  varying  the  number  of  samples  per  signal  collect  with 
the  spinning  radar  emitter. 
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Figure  4.17  reports  the  RMSE  as  the  SNR  at  the  first  signal  collect  is  varied.  The 
RMSE  decreases  as  the  SNR  increases  until  the  frequency,  attitude,  position,  and  velocity 
knowledge  errors  determine  the  performance  of  the  algorithms,  at  which  point  the  RMSE 
begins  to  level  out.  The  AOA/FOA  MLE  geolocation  algorithm  has  the  lowest  RMSE 
throughout. 
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FOA  MLE 
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Figure  4.17:  Sensitivity  analysis  of  varying  the  SNR  at  the  first  signal  collect  with  the 
spinning  radar  emitter. 
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Figure  4.18  reports  the  RMSE  as  the  frequency  knowledge  error  standard  deviation  is 
varied.  The  RMSE  of  the  AOA  LS  and  AOA  MLE  are  constant  since  frequency  estimates 
are  not  included  in  those  algorithms,  while  the  RMSE  increases  for  the  FOA  MLE  and 
AOA/FOA  MLE  algorithms.  The  AOA/FOA  MLE  RMSE  converges  to  the  AOA  MLE 
RMSE  since  the  variances  on  the  angle  estimates  are  lower  than  the  variances  on  the 
frequency  estimates. 
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Figure  4.18:  Sensitivity  analysis  of  varying  the  frequency  knowledge  error  standard 
deviation  with  the  spinning  radar  emitter. 
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Figure  4.19  reports  the  RMSE  as  the  attitude  knowledge  error  standard  deviation  is 
varied.  The  RMSE  of  the  FOA  MLE  is  constant  since  angle  estimates  are  not  included  in 
that  algorithm,  while  the  RMSE  increases  for  the  AOA  LS,  AOA  MLE,  and  AOA/FOA 
algorithms.  The  AOA/FOA  MLE  RMSE  converges  to  FOA  MLE  RMSE  since  the 
variance  on  the  frequency  estimates  are  lower  than  the  variance  on  the  angle  estimates. 
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Figure  4.19:  Sensitivity  analysis  of  varying  the  attitude  knowledge  error  standard  deviation 
with  the  spinning  radar  emitter. 
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Figure  4.20  reports  the  RMSE  as  the  position  knowledge  error  standard  deviation  is 
varied.  The  RMSE  is  relatively  unaffected  for  all  4  algorithms  until  200  m  position  error, 
then  the  RMSE  of  the  FOA  MLE  begins  to  increase.  The  CRLBs  used  do  not  include 
position  errors  and  remain  constant.  The  AOA/FOA  MLE  geolocation  algorithm  has  the 
lowest  RMSE  throughout. 
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Figure  4.20:  Sensitivity  analysis  of  varying  the  position  knowledge  error  standard  deviation 
with  the  spinning  radar  emitter. 
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Figure  4.21  reports  the  RMSE  as  the  velocity  knowledge  error  standard  deviation  is 
varied.  The  RMSE  increases  sharply  for  the  FOA  MLE  and  AOA/FOA  MLE  algorithms. 
The  RMSE  of  the  AOA  LS  and  AOA  MLE  algorithms  is  constant  since  they  do  not  depend 
on  the  velocity  of  the  CubeSat.  The  CRLBs  used  do  not  include  velocity  errors  and  remain 
constant.  In  the  case  of  relatively  large  velocity  knowledge  errors,  the  AOA  LS  and  AOA 
MLE  algorithms  have  lower  RMSE  than  the  FOA  MLE  and  AOA/FOA  MLE  algorithms. 
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Figure  4.21 :  Sensitivity  analysis  of  varying  the  velocity  knowledge  error  standard  deviation 
with  the  spinning  radar  emitter. 
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4.2.2  Burst  Communications  Emitter  Sensitivity  Analysis. 

Simulation  results  of  the  sensitivity  analysis  with  the  burst  communications  emitter  are 
reported  in  Figures  4.22  through  4.29.  The  parameter  sensitivities  observed  with  the  burst 
communications  emitter  are  similar  to  the  spinning  radar  emitter.  The  AOA/FOA  MLE 
geolocation  algorithm  consistently  has  a  lower  RMSE  than  the  other  3  geolocation 
algorithms  with  the  exception  of  high  velocity  knowledge  errors. 
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Figure  4.22:  Orbital  pass  variation  with  the  burst  emitter. 
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Figure  4.23:  Number  of  collects  variation  with  the  burst  emitter. 
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Figure  4.24:  Number  of  samples  per  collect  variation  with  the  burst  emitter. 
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Figure  4.25:  SNR  at  the  first  signal  collect  variation  with  the  burst  emitter. 
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Figure  4.26:  Frequency  knowledge  error  variation  with  the  burst  emitter. 
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Figure  4.27:  Attitude  knowledge  error  variation  with  the  burst  emitter. 


Figure  4.28:  Position  knowledge  error  variation  with  the  burst  emitter. 
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Figure  4.29:  Velocity  knowledge  error  variation  with  the  burst  emitter. 
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4.2.3  Continuous  Communications  Emitter  Sensitivity  Analysis. 

Simulation  results  of  the  sensitivity  analysis  with  the  continuous  communications 
emitter  are  reported  in  Figures  4.30  through  4.37.  The  parameter  sensitivities  observed 
with  the  continuous  communications  emitter  are  similar  to  the  spinning  radar  and  burst 
communications  emitters.  The  AOA/FOA  MLE  geolocation  algorithm  consistently  has  a 
lower  RMSE  than  the  other  3  geolocation  algorithms  with  the  exception  of  high  velocity 
knowledge  errors. 
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Figure  4.30:  Orbital  pass  variation  with  the  continuous  emitter. 
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Figure  4.31:  Number  of  collects  variation  with  the  continuous  emitter. 
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Figure  4.32:  Number  of  samples  per  collect  variation  with  the  continuous  emitter. 
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Figure  4.33:  SNR  at  the  first  signal  collect  variation  with  the  continuous  emitter. 
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Figure  4.34:  Frequency  knowledge  error  variation  with  the  continuous  emitter. 
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Figure  4.35:  Attitude  knowledge  error  variation  with  the  continuous  emitter. 
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Figure  4.36:  Position  knowledge  error  variation  with  the  continuous  emitter. 
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Figure  4.37:  Velocity  knowledge  error  variation  with  the  continuous  emitter. 
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RMSE  (m) 


4.2.4  Geolocation  Accuracy  Comparison  of  Emitter  Types. 

Since  the  AOA/FOA  MLE  geolocation  algorithm  has  been  shown  to  consistently  have 
the  lowest  RMSE  of  the  4  geolocation  algorithms,  the  performance  of  the  AOA/FOA 
MLE  for  the  3  emitter  types  is  compared  in  Figures  4.38  through  4.45. 
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Figure  4.38:  AOA/FOA  MLE  emitter  comparison  for  orbital  pass  variation. 
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jure  4.39:  AOA/FOA  MLE  emitter  comparison  for  the  percent  of  collects  along  pass  7. 
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Figure  4.40:  AOA/FOA  MLE  emitter  comparison  for  samples  per  collect  variation. 
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Figure  4.41:  AOA/FOA  MLE  emitter  comparison  for  SNR  at  first  collect  variation. 


❖ 

Burst  MCS 

— 

Burst  CRFB 

« 

Radar  MCS 

Radar  CRFB 

O 

Continuous  MCS 

— 

Continuous  CRFB 

Figure  4.42:  AOA/FOA  MFE  emitter  comparison  for  frequency  knowledge  error  variation. 
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Figure  4.43:  AOA/FOA  MLE  emitter  comparison  for  attitude  knowledge  error  variation. 
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Figure  4.44:  AOA/FOA  MFE  emitter  comparison  for  position  knowledge  error  variation. 
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Figure  4.45:  AOA/FOA  MLE  emitter  comparison  for  velocity  knowledge  error  variation. 
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V.  Conclusion 


This  thesis  presented  four  geolocation  methods  (LS,  AOA,  FOA,  and  AOA/FOA)  to 
estimate  the  position  of  a  stationary  RF  emitter  from  AOA  and/or  FOA  measurements  at  a 
single  moving  receiver  platform.  A  single  emitter  with  no  co-channel  interference  was 
assumed  to  simplify  the  analysis.  The  MUSIC  algorithm  was  used  for  AOA 
measurements  and  the  frequency  MLE  was  used  for  FOA  measurements.  The  analysis 
scenario  considered  consisted  of  a  single  6U  CubeSat  receiver  platform  in  LEO  receiving 
RF  signals  from  a  terrestrial  emitter.  A  simulation  framework  was  developed  to  validate 
the  statistical  performance  of  the  geolocation  algorithms  against  the  respective  CRLBs 
and  to  conduct  a  system  parameter  sensitivity  analysis. 

5.1  Overall  Research  Conclusions 

From  the  system  parameter  sensitivity  analysis  results  reported  in  Chapter  IV,  the 
AOA/FOA  MLE  geolocation  algorithm  consistently  has  the  lowest  RMSE  of  the  four 
geolocation  algorithms  analyzed  in  this  thesis.  The  increased  performance  of  the 
AOA/FOA  algorithm  is  attributed  to  the  greater  number  of  measurements  available  per 
signal  collect  (angles  and  frequency  estimates)  and  the  intersection  of  the  AOA  and  FOA 
covariance  ellipsoids.  As  observed  in  Figures  3.9  and  3.17,  the  AOA/FOA  covariance 
ellipsoid  is  the  intersection  of  the  AOA  and  FOA  covariance  ellipsoids.  If  an  AOA  payload 
is  present  on  a  single  moving  receiver  platform,  the  implementation  cost  of  incorporating 
FOA  measurements  is  relatively  low  for  an  increase  in  geolocation  accuracy. 

Conducting  the  system  parameter  sensitivity  analysis  in  terms  of  the  SNR  and  number 
of  samples,  rather  than  a  specific  received  signal,  and  parameter  knowledge  error  provides 
intuition  on  the  dependence  of  the  geolocation  algorithms  performance  on  the  various 
system  parameters.  In  general,  the  geolocation  accuracy  increases  as  the  SNR,  number  of 
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collects,  and  number  of  samples  per  collect  increase.  The  geolocation  accuracy  decreases 
as  the  amount  of  frequency,  attitude,  position,  and  velocity  knowledge  error  increases.  At 
high  values  of  SNR  and  number  of  samples  per  collect,  the  geolocation  accuracy  of  all 
four  algorithms  is  determined  by  the  parameter  knowledge  errors.  At  high  values  of 
attitude  and  frequency  knowledge  error,  the  AOA/FOA  algorithm  accuracy  approaches  the 
AOA  or  FOA  accuracy  due  to  using  the  individual  measurement  variances  as  weighting 
factors  in  the  MLE  geolocation  algorithms. 

5.2  Recommendations  for  Future  Work 

Incorporating  multiple  receiver  platforms  into  the  simulation  framework  would  allow 
for  the  analysis  of  additional  geolocation  algorithms  such  as  time  difference  of  arrival, 
frequency  difference  of  arrival,  and  direct  position  determination.  The  geolocation  CRLBs 
and  MLEs  can  be  revised  to  include  the  position  and  velocity  knowledge  errors. 

Additional  analysis  fidelity  can  be  added  by  incorporating  implementation  specific 
considerations  such  as:  mutual  coupling  effects  and  variable  gain  pattern  of  the  receiver 
antenna  array;  variable  gain  pattern  of  the  transmit  antenna;  phase  coherence  and  noise 
characteristics  of  the  RF  receivers;  and  characteristics  of  a  specific  RF  emitter.  Methods 
for  the  data  association  and  segregation  of  multiple  emitters  was  not  addressed  in  this 
thesis  and  the  co-channel  interference  will  decrease  the  performance  of  signal  parameter 
estimates.  Consideration  of  a  moving  emitter  encourages  time  varying,  target  tracking, 
and  motion  analysis.  Incorporation  of  a  surface  of  the  Earth  constraint  with  terrain  data  is 
also  recommended  for  future  research. 
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